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ABSTRACT 


A key drawback to estimating geodetic and geodynamic parameters over time 
based on satellite laser ranging (SLR) observations is the inability to accurately model 
all the forces acting on the satellite. Errors associated with the observations and the 
measurement model can detract from the estimates as well. These “model errors” 
corrupt the solutions obtained from the satellite orbit determination process. Dynamical 
models for satellite motion utilize known geophysical parameters to mathematically 
detail the forces acting on the satellite. However, these parameters, while estimated as 
constants, vary over time. These temporal variations must be accounted for in some 
fashion to maintain meaningful solutions. 

The primary goal of this study is to analyze the feasibility of using a sequential 
process noise filter for estimating geodynamic parameters over time from the Laser 
Geodynamics Satellite (LAGEOS) SLR data. This evaluation is achieved by first 
simulating a sequence of realistic LAGEOS laser ranging observations. These 
observations are generated using models with known temporal variations in several 
geodynamic parameters (along track drag and the J 2 , / 3 , / 4 , and J 5 geopotential 
coefficients). A standard (non-stochastic) filter and a stochastic process noise filter are 
then utilized to estimate the model parameters from the simulated observations. 

The standard non-stochastic filter estimates these parameters as constants over 
consecutive fixed time intervals. Thus, the resulting solutions contain constant 
estimates of parameters that vary in time which limits the temporal resolution and 
accuracy of the solution. The stochastic process noise filter estimates these parameters 
as correlated process noise variables. As a result, the stochastic process noise filter has 
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the potential to estimate the temporal variations more accurately since the constraint of 
estimating the parameters as constants is eliminated. 

A comparison of the temporal resolution of solutions obtained from standard 
sequential filtering methods and process noise sequential filtering methods shows that 
the accuracy is significantly improved using process noise. The results show that the 
positional accuracy of the orbit is improved as well. The temporal resolution of the 
resulting solutions are detailed, and conclusions drawn about the results. Benefits and 
drawbacks of using process noise filtering in this type of scenario are also identified. 
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CHAPTER 1 


INTRODUCTION 


1.1 Overview 

The current state of the art in filtering Earth orbiting satellite data has reached the 
point where temporal variations in the gravity field (particularly J 2 and 7 3 ) appear to be 
observable. Determining these variations is of interest for determining global changes 
in mass distribution as well as for insight into interior mass properties. The desire to 
then obtain accurate estimates of these variations, as well as temporal variations in other 
geophysical parameters, provides the motivation for this study. In particular, it is of 
interest to determine if the relatively sparse, but accurate, laser range tracking of the 
LAGEOS satellite can be used to resolve these variations in the low degree coefficients 
of the Earth’s gravity field with the use of a stochastic filter. Conventional filtering 
methods typically estimate these types of variations as piecewise constants over a given 
data arc [Nerem et al., 1993]. This results in a discontinuous solution for the variations 
with limited temporal resolution (e.g. monthly). Stochastic filtering methods can 
estimate these variations as continuous process noise parameters which are correlated in 
time. This type of estimation procedure has the potential for resolving these variations 
much more accurately. This study analyzes the ability of a stochastic filter to recover 
and accurately estimate such variations. 

1.2 Model Errors 

A dynamical model is a mathematical representation of the forces acting on a 
physical system. If the dynamics of the system are known perfectly, then a dynamical 
model for the system can be used to determine the exact physical state of the system for 
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all time. Most physical systems cannot be modeled or observed perfectly, as is the case 
with Earth orbiting satellites. Errors in the dynamical model for a satellite result in 
differences between the true satellite state and the predicted satellite state based on the 
dynamical model. Many of these dynamical model errors are a result of not accurately 
knowing the values of physical parameters which help define the satellite dynamics 
mathematically. Geopotential coefficients, drag coefficients, and solar radiation 
pressure coefficients are some examples. Also, model errors may consist of 
unpredictable time variations in the model parameters which can not be modeled. 

Likewise, errors may exist between the actual observation or measurement and 
the model used to determine a computed or predicted measurement based on the satellite 
state. These are measurement model errors. A computed measurement of the state is 
necessary to compare to the actual observed measurement of the state in order to 
compare the current actual state with the current predicted state. Errors due to moving 
laser stations (plate tectonics), clock errors, and atmospheric refraction are some 
examples of measurement model errors. There are also some model errors, such as 
polar motion errors, which may be considered both dynamical model errors and 
measurement model errors as these parameters are present in both models. 

1.3 Orbit Determination Filters 

By observing the satellite over time, the differences between the true (observed) 
state and the predicted state based on the dynamical model become apparent. Orbit 
determination is the process of obtaining the best estimate of the state of a satellite based 
on observations of that satellite. An orbit determination filter is used to combine the 
information from a set of observations into an estimate of the state while filtering out 
the errors associated with the observations and the dynamical model. 
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There are several types of orbit determination filters. The particular application 
usually dictates the specific filter which is most appropriate. For this study, the two 
main types of filters which are of concern are a standard (non-stochastic) Square Root 
Information Filter (SRIF) and a process noise (stochastic) SRIF. These two filters are 
compared and their ability to resolve specific model errors is assessed. 

The satellite state parameters may be defined as the position and velocity of the 
satellite plus constant model parameters. These model parameters are defined as 
constants in the mathematical model, yet they may in fact vary in time. Such 
parameters must be defined as constants in the model since any variations are usually 
unpredictable. If variations of a particular model parameter are predictable, then these 
variations can be built into the model. The resulting dilemma is that of estimating 
unpredictable variations in model parameters as constants. The limitations inherent to 
this approach are obvious. 

Standard non-stochastic filters estimate model parameters (used in the 
dynamical and measurement models) as constants. Process noise stochastic filters can 
estimate these model parameters as stochastic, time varying parameters correlated in 
time. Generally, stochastic filters are most often used to estimate extra accelerations 
which account for the total effect of the individual model errors for various parameters 
combined [Yunck et al., 1990; Wu et al., 1992; Yunck et al., 1994]. These stochastic 
acceleration estimates result in an estimated state which is much closer to the true state 
(position and velocity). The estimates of the constant parameters of the state that are 
used in the dynamical and measurement models are generally no closer to the true 
values, however. This is due to the fact that the errors in these parameters have not 
been estimated, but rather the net effect of these errors has been estimated as additional 
accelerations. This study focuses on the ability of a stochastic orbit determination filter 
to estimate time varying model parameters directly and simultaneously with the satellite 
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state, and the associated accuracies of the estimates of the model parameters and the 
satellite orbit. 

1.4 Satellite Laser Ranging and LAGEOS 

Satellite laser ranging (SLR) is a satellite observational technique which has 
developed and matured over the last couple of decades. The observable used in SLR is 
a range measurement from a ground based laser station to the satellite. This range 
measurement is based on the round trip travel time of a laser pulse and the constancy of 
the speed of light. Current precisions for SLR measurements are at the sub-centimeter 
level [. Kolenkiewicz et al., 1991]. This highly accurate type of measurement results in 
geodynamic solutions which are more accurate than those solutions using other types of 
ground-based measurements. 

Launched in 1976, the Laser Geodynamics Satellite (LAGEOS) has been a 
popular target for SLR. Due to its high altitude (approximately 5900 km), the 
LAGEOS orbit is highly predictable and fairly easy to model [Cohen and Smith, 1985]. 
The high altitude results in smaller effects from atmospheric drag or the higher 
frequency geopotential terms. Thus, LAGEOS solutions are less sensitive to model 
errors from drag and gravity than lower altitude satellites, and more accurate 
geodynamic solutions are possible. 

1.5 Thesis Objectives 

The main objective of this thesis is to assess the ability of a sequential process 
noise filter to resolve specific model deviation signals embedded in simulated LAGEOS 
SLR data. The accuracies of the estimated model deviation signals and the satellite state 
are analyzed by comparing them to the known truth used in the simulated orbit. 
Specifically, model deviation signals in an along track drag parameter and the second 
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and third degree zonal geopotential coefficients (7 2 and 7 3 ) are introduced into a 
simulated one year orbit and temporally resolved by a stochastic filter. Comparisons 
are made to a solution generated with a standard non-stochastic filter. Also, solutions 
are generated in the presence of additional model deviation signals (7 4 and 7 5 ), which 
are not estimated. A simulated three year orbit with model deviation signals present in 
C„ 7 2 , and 7 3 is also processed just using the stochastic filter. The three year arc is not 
processed using the standard non-stochastic filter. The ability of a stochastic process 
noise filter to resolve variations over this longer three year arc is assessed. The 
advantages and disadvantages of using a stochastic filter to generate a LAGEOS 
geodynamic solution are addressed as well. 

This research will benefit future geodetic and geodynamic studies which attempt 
to temporally resolve satellite model parameters. The feasibility of this approach is 
assessed as it is related to orbit determination and satellite geodynamics. The specific 
benefits and drawbacks of sequential process noise filtering in determining these 
solutions will help to refine the role of this type of filtering in future analyses. This 
study will also benefit orbit determination studies considering the use of process noise 
parameters and stochastics in general in estimating model errors. 

1.6 Description of Chapters 

Chapter two reviews the basic aspects of filtering theory. The orbit 
determination problem is presented and state estimation theory is summarized. The 
square root information filter (SRIF) is developed both with and without correlated 
process noise. The standard SRIF (without process noise) and the stochastic (process 
noise) SRIF are the two orbit determination filters compared in this analysis. 

Chapter three discusses the satellite laser ranging technique (SLR) and the Laser 
Geodynamics Satellite (LAGEOS) and their roles in geodynamics. Some current 
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LAGEOS geodynamic solutions are presented and possible benefits of stochastic 
filtering of LAGEOS data are detailed. 

Chapter four details the simulation model used in the generation of the simulated 
LAGEOS data. The measurement model and tracking station network are presented 
along with the dynamical model. The specific model deviation signals inserted into the 
simulated data are presented as well. 

Chapter five details the results obtained from filtering the simulated data with 
both the standard and the stochastic filter. The temporal resolution of the model 
deviation signals estimated by each filter are compared. The accuracies of each filter in 
estimating the satellite orbit are also compared. 

Chapter six discusses the conclusions drawn based on the results of filtering the 
simulated data. Benefits and drawbacks of using process noise to estimate variations of 
parameters from their nominal values are discussed. Finally, recommendations for 
future studies related to this research are made. 



CHAPTER 2 


FILTERING THEORY 


2.1 Introduction 

By observing a physical system, an estimate of the state of the system can be 
made based on the observations. The state of the system defines what the system is 
doing or looks like at any given time. Generally, the measurements or observations of 
the system usually contain some type of noise. Filtering or estimation is the process of 
detennining the best estimate of the state of the system (by some measure) from these 
noisy observations. By “filtering” out the noise on the observations, a “best estimate” 
of the state can be made which is often more accurate than the noise on the 
observations. 

In the process of detennining the best estimate of a satellite state, observations 
of the satellite are required. These observations are processed in such a way as to filter 
out the errors associated with them as the state of the satellite is estimated. This chapter 
describes the theory associated with this process, which is commonly referred to as 
filtering. 

2.2 The Orbit Determination Problem 

If the forces acting on a satellite have been modeled perfectly and the initial 
conditions of the satellite are known exactly, then the state of the satellite can be found 
for all time by integrating its equations of motion. A general form for the equations of 
motion is 


r = F(r,r,p,0 


( 2 . 1 ) 
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where r and r are the geocentric position and velocity of the satellite, p are constant 
parameters contained in the mathematical force model for the satellite’s dynamics, and t 
is the integration time. In general, however, neither the dynamical force model nor the 
true initial conditions are known exactly. Perfect dynamical force models do not exist, 
and the constant parameters p in the force model have errors as well as possibly 
unpredictable changes over time. True initial conditions usually differ from those 
calculated a priori, so even a “perfect” dynamical force model for the satellite will 
predict a trajectory that differs from the true trajectory as a result of this initial condition 
error. Also, some accelerations may not be modeled at all by the chosen dynamical 
force model. As a result, observations during the satellite’s orbit must be made in order 
to determine or verify its subsequent trajectory. Since the state variables associated 
with a satellite ( r and r ) cannot be observed directly, other measurements such as a 
range or a range rate must be made. The satellite state can then be estimated based on 
these measurements. This is known as the orbit determination problem. 

Due to inherent measurement errors in making observations of the satellite, the 
trajectory estimated based on the observations will be different than the satellite’s true 
trajectory, even if the dynamical force model for the motion of the satellite is perfect. 
Errors in the measurement model, which relates the satellite observations to the satellite 
state, also exist. The observation vector Z of the satellite at time t is usually a nonlinear 
function of the satellite state ( r and r ), a set of measurement model parameters b, plus 
some random measurement noise e z : 

Z = G(r,r,b,0 + e z (2.2) 

Thus, the satellite state and the observations are related in a nonlinear manner, denoted 
by the function G. By linearizing this nonlinear function about a known reference 
trajectory, the orbit determination problem can be simplified [Smith et al., 1962]. If the 
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reference trajectory and the true trajectory are sufficiently close during the time interval 
of concern, then the actual trajectory can be expanded in a Taylor series about the 
reference trajectory for each point in time. By truncating this expansion after the first 
order terms, the state deviation from the reference trajectory may be represented by a set 
of linear differential equations. Likewise, the observation deviation may be linearly 
related to the state deviation. Thus, the nonlinear problem of determining the satellite’s 
state can be transformed to a linear problem which determines the satellite’s state 
deviation from some reference state (trajectory). 

To summarize mathematically, at time t, the deviations in the satellite state and 
observations are: 

x(t) = X(t)-X\t) (2.3) 

z(0 = Z(r) - Z\t) (2.4) 

where X(f) is the true satellite state; X*(f) is the reference state based on a specified 
dynamical force model; x(r) is the state deviation from the reference state X’(r); Z (t) is 
the observation vector; Z *(r) is the computed observation vector based on the reference 
state X*(f); and z(t) is the observation deviation from the computed observation Z *(r). 
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By using the deviations defined in equations 2.3 and 2.4, and dropping the time 
dependence of the state and observation deviation for notational simplicity, equations 
2.1 and 2.2 may be linearized as 


x = A(t)\ 


(2.5) 


z = Hx + e z 


( 2 . 6 ) 


where 


3F . 

A(t) = ^ X,/) 


(2.7) 


ff=ff<xV) 


(2.8) 


For the nx 1 state vector X, an nxn state transition matrix <J> is defined 




dX(t) 

dX(0 


with 


0(t,t / ) = A(t)O(t,t / ) , 0(^,0 = I 

and 

Thus, the linearized system (equation 2.5) has the general solution 


(2.9a) 


(2.9b) 


(2.9c) 


x(0 = 0(t,tj)x(tj ) 


( 2 . 10 ) 
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If desired, by using this result and choosing an arbitrary epoch time t a , the linearized 
observation-state relation (equation 2.6) may also be written 

z = H 0 X 0 +e z (2.11) 

where 

H 0 = H®(t,t 0 ) and x 0 = x(/ 0 ) (2.12) 

For anmxl observation deviation vector z, Hisamxn matrix. Thus, observation 
deviations at any time t are linearly related to state deviations at that time or at some 
epoch time t a . 

2.3 State Estimation 

Least squares filters are the most commonly used filters for orbit determination 
problems since they provide the best estimates of the state when the uncertainty is due 
to Gaussian or normally distributed noise. They minimize the square of the difference 
between the observed measurement and the expected measurement computed from an 
observational model based on the state of the reference trajectory. 

Observations can be filtered in two basic modes: batch or sequential. 

Depending on the specific problem and application, data is filtered as an entire batch, 
sequentially, or a combination of both. A batch filter estimates the state at an epoch 
time based on all the observations taken over a given time interval. A sequential filter 
processes the observations one by one, estimating either the current state or some epoch 
state after each observation. By processing the observations sequentially, many large 
matrix inversions and multiplications are avoided since the size of the observation 
deviation vector z is dependent only on the number of observations at the current time 
rather than the number of observations over the entire time interval. When computer 
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storage is limited, the batch method is less favorable due to the necessity to store all the 
observations and iterate through the data until a solution is converged upon. Batch 
methods, however, are easy to implement and less sensitive to erroneous data points. 
While still quite simple to implement, sequential filters are more sensitive to erroneous 
data and numerically more unstable. Sequential algorithms are often used in 
applications where computer storage is a limitation. 

In general, a “standard” filter is most frequently used to estimate the state of 
satellites and model parameters. Tapley et al. [1993] and Nerem et al. [1994] are two 
typical examples of recent satellite solutions that estimate geodynamic parameters. 
Nearly all current satellite solution methods use some type of standard least squares 
orbit determination filter to estimate satellite states and model parameters. When epoch 
values for the satellite position and velocity are used, a standard filter estimates all of 
the parameters in the state vector as constants, which are inherently assumed to be time 
invariant. The epoch satellite position and velocity define the satellite state at a specific 
epoch time, and the model parameters may be defined in either the dynamical force 
model (such as geopotential coefficients or drag coefficients) or the measurement model 
(such as tracking station positions). Errors inevitably enter into the solution or estimate 
when model parameters are constants in the dynamical or measurement model, but in 
reality, vary in time by some measurable, yet unpredictable amount. This type of error 
is often referred to as model error. That is, the system is not performing or behaving 
like the predicted model. This problem can be minimized somewhat by using multiple 
piecewise constants which represent the particular variation as a group of consecutive 
discontinuous constants, each representing a shorter time interval than the entire time 
period of interest. The effectiveness of this approach is dictated by many properties, 
including the observability of the variation the type, amount, and density of available 
data, and the frequency of the variation. 
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Another type of filter, called a “process noise” or “stochastic” filter, can be used 
to address this type of error. A process noise filter estimates specified “constant” 
parameters as stochastic, time varying parameters. Lichten [1990a, 1990b] details 
procedures that involve process noise for estimating geodynamical parameters using 
Global Positioning System (GPS) data. Lichten [1992] and others have presented 
solutions which involve the use of process noise parameters in satellite orbit estimation. 
Even though the parameters are modeled as constants in the dynamical force model or 
the measurement model, the filter allows the estimates for these parameters to vary in 
time. The parameters can vary in a manner which is consistent with the information 
contained in the observations. Examples of phenomena which can be approximated 
quite well by process noise parameters are solar radiation pressure, mismodeled drag 
effects, leaky attitude control systems, moving tracking stations, polar motion and 
Earth rotation parameters, clock errors, atmospheric path delays, gravity field model 
errors, and linearization errors. Thus, a process noise filter is capable of reducing or 
eliminating many types of model error. 

Often, process noise filters are utilized to lump all model errors into additional 
acceleration parameters or a specific model parameter to reduce the negative effect of the 
total model error on the estimate of the orbit. Reduced dynamic tracking techniques 
have utilized this approach [Yunck et al., 1990; Wu et al., 1991; Yunck et al., 1994; 
Gold, 1994]. This approach is used when the temporal variation of model parameters, 
or perhaps their actual values, is of little interest and the objective is only to improve the 
continuous estimate of the satellite position and velocity over time. By representing the 
cumulative effect of all model errors with an additional stochastic acceleration, this can 
be accomplished, and is referred to as the reduced-dynamic •technique [Yunck et al., 
1990]. 



In contrast, estimates of particular parameters may be made stochastically to 
understand their temporal variations in addition to reducing the effect of the model error 
on the estimate of the orbit. This study focuses on the feasibility of this approach in 
accurately estimating some basic geodynamic parameters from LAGEOS data. Thus, 
eliminating various model errors associated with specific parameters through process 
noise filtering may lead to more accurate estimates of the satellite orbit, as well as better 
temporal estimates of the particular dynamical model parameters and measurement 
model parameters. 

2.4 Standard Filtering 

The determination of the best estimate of x satisfying the linearized observation- 
state relation (equation 2.6) is discussed. The least squares best estimate of x is that 
value which minimizes the sum of the squares of the computed observation residuals 
e z . A performance index J is defined 

J = e z e z = II e z II 2 (2.13) 

which must be minimized in a least squares sense to determine the best estimate of x . 
The notation II • II is used to express the Euclidean norm of a vector ( II a 11= V(a T a) ). 
Based on the observations z and specified values of x, the squares of the observation 
errors e z may be summed and minimized. Using equation 2.6, the performance index J 
may be written as 

J = (z - Hx) 7 (z - Hx) (2.14) 

This performance index is now minimized with respect to the state deviation x: 
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or 

H t Hx = H t z (2.16) 

Equation 2.16 is often referred to as the normal equations. The value of x which 
minimizes the performance index J, and thus satisfies the normal equations, is the least 
squares best estimate of the state and denoted x . The error covariance matrix P that is 
associated with the best estimate x is defined 

P = (H t H)~' (2.17) 

so 

x = PH x z (2.18) 

The error covariance matrix P is updated to the time of the new observations, if 
necessary (for current state filters), as follows: 

P,=^(hA)P,^\t2A) (2-19) 


2.5 Standard SRIF Filtering 

To achieve more stability, accuracy, and better numerical conditioning, 
factorized sequential filtering algorithms can be used [ Bierman , 1977]. By updating a 
factorized variation of the covariance matrix for the system, the algorithm exhibits 
improved numeric behavior with less sensitivity to divergence or erroneous data points. 
This is mainly a result of reducing the numerical ranges of the values of the covariance 
matrix (10 -n to 10" => 10 _n/2 to 10 n/z ). Square root factorizations, first introduced by 
Potter (cf. Battin, 1964), are common and quite useful in sequential filter algorithms. 
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A current state sequential square root filter is used in this study, and is described 
below. 

The basic aspects of a commonly used square root filter, known as the square 
root information filter (SRIF), are outlined. In a SRIF, the square root of the 
information matrix associated with the linearized system is operated on rather than the 
covariance matrix. The complete formulation for the SRIF is detailed in Lawson and 
Hanson [1974] and Bierman [1977]. 

If the linearized system deviation defined by equation 2.10 has an associated 
error covariance matrix P, then its associated information matrix A, is defined as 

A = r' = R t R (2.20) 

where R is the square root of the information matrix and is upper triangular in form. 
The a priori error covariance and estimate [p x] are related to the a priori information 
array [ R z] as follows: 

Since 

P = {R T RY' = R-'R~ r (2.21) 

the normal equations (equation 2. 16) may be written 

x ={R T R)~ ] R T z =/T'z (2.22) 

Thus, the best estimate x must satisfy 


Rx = z 


(2.23) 
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So R and z , which represent a priori values of R and z, correspond to the following 
data equation for the system 

z = Rx + e z (2.24) 

where R is square, x is the state of the system, and e z has zero mean with unit 
covariance. The least squares solution to this a priori data equation and the equation for 
the new observations (equation 2.6) is desired. 

The matrix R is updated to the time of the new observations, if necessary (for current 
state filters), as follows: 


R h = (2-25) 

By applying the orthogonal transformation method [Givens, 1959; Householder, 1964; 
Schmidt, 1967] to equations 2.24 and 2.6, the least squares solution of 


R 

H 




(2.26) 


is also the least squares solution of 


T 


R 

H 


x = T 


z 

z 



(2.27) 


where T represents the orthogonal transformation matrix which preserves unit 
covariance characteristics of the error terms. The orthogonal transformation matrix T 
can be chosen such that 


R 


' /s " 

R 

H . 


.0 _ 


(2.28) 


with R upper triangular. 
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Using the same T, the following can be defined: 


z 


■ /V ' 

z 

and T 

e 2 


/s 

e z 

. z . 


- e ls. 


. e 7 


. e Cb _ 


(2.29) 

As defined, e, s represents the error in the least squares fit. Now equation 2.26 becomes 

(2.30) 


R 

0 


x = 


z - e z 

L e Is ~ e e 


Thus, the augmented observation array 


R z 
H z 


(2.31) 


is transformed to 


R z 

0 e, s 


(2.32) 


and observations can be accumulated sequentially with the a priori information array. 
After the final observation is accumulated (transformed), the least squares best estimate 
of the state x , is: 


/N __ 

x = R z 


(2.33) 


This estimate is easily computed via backsubstitution rather than by matrix inversion 
since R is upper triangular. The filtered best estimate of the state deviations x , can be 
computed at any given time. The filtered best estimate of the state is thus based on all 
observations up to that particular time. Smoothed estimates can be computed for all 
times by mapping the state deviations to any previous time. Alternatively, the estimated 
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state deviations can be added to the reference state, and the estimated state at other times 
can be found by integrating the equations of state. 

2.6 Process Noise SRIF Filtering 

By using a process noise filter, many phenomena, unmodeled and mismodeled 
by the dynamical model or measurement model, can be accounted for in a more 
reasonable way. These unmodeled and mismodeled phenomena degrade the solution 
accuracy and are reflected in the measurement residuals (the difference between the 
measured observation and the estimated observation resulting from the solution). 
Process noise filtering can help prevent these model errors from corrupting the 
solution. The development of this sequential current state algorithm is based on the 
SRIF algorithm with correlated process noise, first introduced by Andrews [1968], and 
presented by Bierman [1977]. 

First order exponentially correlated process noise is often called colored noise. 
The process noise filter developed is based on this type of process noise, but other 
types of process noise can be incorporated into process noise filters. Colored noise can 
be described mathematically as: 

± = -(ljp + a (2.34) 

where p is the process noise parameter, T is the time constant of the process, and (O is 
white noise with zero mean. Equation 2.34 can be converted to a recursive form, 
which is useful in discrete time systems, such as the orbit determination problem: 


Pj+l = m jPj + Wj 


(2.35) 
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with 


mj = exp[— (/*, - /,)/ t] (2.36) 

and 


exp[-(^, - (2-37) 


It can be shown that the variance q, associated with the process noise w, is: 


q = {\-m])(? 


(2.38) 


where o 2 is the steady state variance associated with p. 

The state is now partitioned into process noise, dynamic, and bias parameters 
denoted by p, x, and y respectively. The bias parameters are simply constant (time 
invariant). The dynamical model for this partitioned, linearized state is: 


p 

X 


M 0 0 

% % 

P 

X 

+ 

' w / 

0 

y 

M 

0 0 / 

_y . 

/ 

0 


This corresponds to equation 2.10 for the standard SRIF development. The state 
transition matrices & p , <t> x , and & y map state elements at time ( to time . / is a square 

unit matrix equal in size to the number of bias parameters y. M is the process noise 
transition matrix and is assumed diagonal, with diagonals m, given by equation 2.36. 
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The value of m is dictated by the time constant for the process T, in the following way: 


m = 0 T— > 0 [white noise process - not correlated in time] 

0<m<l 0 < T< 1 [colored noise process - somewhat correlated in time] 

m = 1 r — » °° [random walk process - strongly correlated in time] 

The behavior of the process noise parameters p are controlled mainly by the mapping 
elements m associated with each particular parameter. Each of the process noise 
parameters p can have an independent time constant r (and corresponding m) and 
steady state standard deviation o associated with it. A parameter behaves as white 
noise when m = 0, and as a random walk when m = 1 . A parameter behaves like 
colored noise when m is between these two extremes. When a parameter is defined to 
model a random walk process (m = 1), the variance q for that particular parameter must 
be explicitly defined since no steady state exists and the steady state variance & for that 
parameter is undefined. Figure 2.1 shows a white noise and a random walk process 
generated using equation 2.35. Both processes shown in Figure 2. 1 use an initial value 
of p - 0 for the process noise parameter and the same process noise series w (zero mean, 
unit covariance) with a specific value of m. 
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Paralleling the development of the standard SRIF, the mathematical 
development of determining the least squares best estimate using the process noise 
SRIF algorithm follows (this development follows that of Bierman [1977]): 

At time tj, the data equation for R and z can be written: 


/N A. /S 



A A A 



0 0 R y 




(2.40) 


R v is zero unless no data was processed at the last time update. If R is not at the time 
of the current observation (f y+1 ), it must updated. The time update for R is a bit more 
complex in the process noise SRIF algorithm than that described by equation 2.25 in 
the standard SRIF algorithm. From equation 2.39: 


x, +l = & x Xj + O p p y + <P y y 


(2.41) 
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Solving for x, in terms of x ;+1 and substituting this into equation (2.40): 

R p - (R px ® x )® p R px <t>:' R py - (R px & x )<t> y I - P; 1 \z p 

R xp - (R x ® x )& p R x 0 x R^-CR^M, Xj+i = z x (2.42) 

o o R y JL y J [%_ 

But, p ;+1 is needed in equation 2.42 as well. By using the a priori data equation for w, 
this can be done: 

R„v?j = z w (2.43) 

where 

vq 

with R w diagonal and q is the variance associated with the process noise w y , as 
defined in equation 2.38. 

And from equation 2.39 

W; = P; + , - Mpj (2.44) 

so 

RJVw ~ M Pj) = (2.45) 

Now combining equation 2.45 with equation 2.42 to obtain the updated R and z at the 
next time, (r+l : 

-RJ4 R w 0 0 1[ Py 1 \zj 

P /+1 /^. 

(R p -R px %) 0 R px (R py -R px <t> y ) Xj+i = z p (2.46) 

(R xp -R x <t> p ) 0 R x (R xy -R x <t> y ) JL y J [z x 
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where 

R px = R P &' and R x = R/b x (2-47) 

Now, using an orthogonal transformation to upper triangularize [/? z] and eliminate p/. 

-R„M R w 0 0 z w R' P R' pp R' px R'py z p 

T (R p - R px <Z> p ) 0 R px (R py -R px <t> y ) z p = 0 R„ R px R py z p (2.48) 
(Rtp - R x <t> p ) 0 R x (R^-Rf D y ) zj [0 i? xp £ x R xy z x 

Since y is a bias, it is unaffected by the mapping from /, to O+i and need not be included 
in the time update for R and z or in the orthogonal transformation: 

[r,z, ]; = [£,?,];*, (2.49) 

The quantities superscripted by * are not used to obtain filtered estimates, but are critical 
if smoothing is to be done following the filtering. The ~ quantities from the previous 
time are now represented by the ~ quantities at the current time, completing the time 
update for R and z. 
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Next, the measurement accumulation is written: 


Rp Rpx Rpy 

R xp Rx R X y 

0 0 /?, Zy 


T 


H 


z 


Rp Rpx Rpy z P 

0 Rx Rxy Z, 

0 0 Ry Zy 

0 e, s 


(2.50) 


where z is the observation deviation vector for the current time, H is the observation 
partials matrix, the R and z partitions represent the updated R and z at the current time 
after accumulating the current observation, and e ls represents the sum of squares of the 
error terms. At any given time, the best estimate is: 

x = R‘z (2.51) 


where 



r d " 
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Again, since R is kept in upper triangular form, a matrix inversion is not 
necessary to compute the best estimate of the state deviations. The best estimate, x , 
can be found directly via backsubstitution. Similarly, the error covariance matrix P 
for the estimate x is 

P = R R (2.52) 

This best estimate of the state deviations x , is referred to as a filtered estimate. 
It represents the best estimate at that particular time based on all observations up to and 
including the observation at that time. Smoothed estimates represent the best estimate at 
a given time based on all of the observations, both before and after the time of the 
estimate. Smoothing is the process of determining these smoothed estimates once all 
the observations have been processed. 

If these smoothed estimates of the state at other times are desired, then the 
smoothing parameters saved at each time (from equation 2.48) may be used. Given 

[p ; , Xp y ; ], then [p ; _,, x ; _, , y ; _,] can be found as follows. Since y is a constant, its 

smoothed estimate equals its filtered estimate: 

yj-i = y, (2.53) 

Using equations 2.46 and 2.48, the smoothed estimate of p is 


_ _ Kl - pAKpI - x A R lx\j - y AKyl 

ia 


(2.54) 


and using equation 2.39, the smoothed state estimate is 

Xy_, = ST' [xj - <D p py_, - O y y hl ] 


(2.55) 
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Using these recursive equations, smoothed estimates of the state can be found for all 
previous times. This smoothing process is slightly more complex when stochastic 
parameters are present than smoothing for the standard SRIF which only involves 
integrating the equations of state backwards or mapping state deviations using the state 
transition matrix <I>. 

Filtered estimates of the error covariance matrix P at any given time are found 
by using equation 2.52. Once all the observations have been processed, the error 
covariance matrix at previous times may be desired (smoothed error covariance values). 
This requires the knowledge of either the information matrix R at any previous time so 
that P may be found via equation 2.52 or the knowledge of P directly by way of 
mapping from the final time. The mapping of the information matrix R, based on R at 
the final time after all observations have been processed, may be done by using the time 
update for R backwards. While equations 2.46 and 2.48 are used to update R forward 
in time during the filtering, they may be used for backwards mapping by substituting 
t H for t j+ 1 in calculating the state transition matrices. Given the final information matrix 
R at the final time based on all observations (equation 2.51), equations 2.46 and 2.48 
may be used to propagate R to previous times. Then, the ~ quantities resulting from 
equation 2.48 represent R at the previous time rather than the current time, and 
the ~ quantities represent the current time rather than the previous time. 

Alternatively, P from the final time may be mapped directly to other times by 
using the generalized error covariance update relation: 
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Recall, the state transition matrices & p , d> x , and <&>. map state elements at time r, to time 
M is the process noise transition matrix and is assumed diagonal with diagonals m 
given by equation 2.36, and Q is the diagonal error covariance matrix for the process 
noise w with diagonals given by equation 2.38. I is a square unit matrix equal in size 
to the number of bias parameters y and I p is a square unit matrix equal in size to the 
number of process noise parameters p. In order to map P to previous times using 
equation 2.56, the substitution of t H for must be made in calculating the state 
transition matrices and Q. 

The process noise SRIF algorithm parallels the standard SRIF algorithm with 
the exception of the propagation of R and the smoothing portion which are both 
different from the standard SRIF algorithm since stochastic parameters (p) are 
involved. The time variation of each process noise parameter is controlled with a time 
constant T, and a steady state variance a 2 , for the parameter. As T increases and a 
decreases, the corresponding parameter stays more constant from one time to the next. 
As T decreases and a increases, the corresponding parameter becomes more time 
varying. These two parameters are “tuned” or chosen depending on the expected time 
variation of the process the process noise parameter is modeling. 

2.7 Conclusions 

The basic attributes of the standard and stochastic filters are now summarized. 
As mentioned previously, the inherent disadvantage of using a standard filter to 
estimate temporal variations of parameters is the fact that a standard filter must estimate 
these variations as constants. In order to minimize this inadequacy, a parameter may be 
divided into multiple “sub-parameters”, each of which represents the original parameter 
for a particular and unique time interval during the period of interest. For example, 
over a one year data arc, the solar radiation pressure coefficient of a satellite may be 
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estimated using 52 consecutive parameters, each representing a unique one week period 
during the data arc. The consecutive joining of these 52 discontinuous sub-parameters 
over the one year time period represents the year long estimate for the parameter. 
Based on the appearance of the consecutive joining of these discontinuous sub- 
parameters, this type of solution is denoted a “boxcar” solution for the purposes of this 
study. For parameters which are not purely kinematic in nature, this involves adding 
each of these sub-parameter boxcars to the satellite state vector. Thus, the 
computational price for filtering a two year arc instead of a one year arc doesn’t merely 
double, but grows geometrically depending on the number of parameters that are 
estimated with shorter boxcars. For data arcs spanning long time periods such as six 
months or more, estimating the variations of parameters in this fashion can be quite 
tedious and computationally laborious. Further, splitting a parameter into multiple 
boxcars does not guarantee improved recovery or resolution of the variation. The ideal 
minimum time interval for a parameter depends on many things, including the 
sensitivity of the satellite to variations in the parameter, data type and density, and the 
observability of the variation in the satellite measurement. As an example, using 
boxcars that span time intervals less than one month for variations in gravity 
coefficients generally does not improve the recovery of these variations. 

In contrast, by estimating variations in a parameter as stochastic process noise 
parameters, no increase in the size of the satellite state vector is needed. Since the 
parameter is inherently stochastic as opposed to constant, no additional sub-parameters 
are required. As shown in the development of the process noise SRIF, stochastic 
estimates of the process noise parameters are obtained at every observation time. If the 
variation in a parameter is recoverable via this approach, the estimation of the variation 
may be significantly simplified computationally. 



CHAPTER 3 


SATELLITE LASER RANGING AND LAGEOS 

3.1 Introduction 

Some basic principles of the satellite laser ranging technique are discussed. The 
specific role of the Laser Geodynamics Satellite (LAGEOS) in satellite laser ranging is 
outlined as well. Some current solutions based on LAGEOS data are also presented. 
Benefits of processing LAGEOS data using process noise filtering techniques are also 
proposed. 

3.2 Satellite Laser Ranging 

Satellite Laser Ranging (SLR) is a highly accurate and precise tracking method 
used for orbit determination of satellites. Ground based tracking stations measure the 
distance (range) from the satellite to the station by using lasers. The time it takes the 
photon pulse from the ground based laser to travel to the satellite and back is converted 
to a range measurement by using the constancy of the speed of this light pulse. In the 
case of satellites designed for laser ranging targeting, the laser pulse is reflected by a 
cluster of cube-comer reflectors located on the satellite surface. The precision of SLR 
is now at the sub-centimeter level [Kolenkiewicz et al., 1991; Tapley et al., 1993; 
Degnan, 1994], and through the formation of laser normal points, can reach a few 
millimeters [ Degnan , 1994] (see Smith et al. [1991] for a detailed discussion on normal 
point formation and processing). The accuracy of SLR ranging systems is also at the 
sub-centimeter level [Degnan, 1985, 1993, 1994], and the number of tracking stations 
around the world has been increasing consistently over the past two decades. By 
integrating these highly accurate and precise SLR measurements from multiple tracking 



31 


stations with a dynamical model for the motion of the satellite, numerous model 
parameters can be determined to a high degree of precision via a least squares orbit 
determination filter. These parameters include those describing the satellite’s motion, 
positions of the tracking stations, tectonic plate motions, both spatial and temporal 
variations of the Earth’s gravity field, Earth rotation and polar motion parameters, Earth 
and ocean tides, variations in the center of mass of the total Earth system, and other 
geodetic parameters used in the dynamic and measurement models [Degnan, 1994]. 
The accuracies of geodetic and geophysical results obtained through satellite 
observations are directly related to, and often limited by, the accuracies of the 
observables used in obtaining the results. The accuracy of the SLR observable has 
enabled the quality of geophysical and geodetic results to improve significantly through 
the use of laser target satellites. 
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3.3 The Laser Geodynamics Satellite 

The Laser Geodynamics Satellite (LAGEOS), launched by NASA from the 
Western Test Range in California on May 4, 1976, is a passive, spherical, artificial 
satellite. Figure 3.1 shows LAGEOS which resembles a large cannonball covered with 
retroreflectors. 



Figure 3.1. The Laser Geodynamics Satellite (LAGEOS), from Smith et al. [1991]. 


As outlined by Fitzmaurice et al. [1977] and Cohen and Smith [1985], Table 3. 1 shows 
the orbit and satellite characteristics of LAGEOS. It is covered with cube-comer 
reflectors, 422 made of fused silica and four made of germanium. It was designed 
specifically and exclusively as a long-term laser ranging target to improve geophysical 
and geodetic solutions through the use of SLR data. 
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Table 3.1. Orbit and Satellite Characteristics of LAGEOS. 


Diameter 

60 cm 

Mass 

411 kg 

Fused silica retroreflectors 

422 

Germanium retroreflectors 

4 

Inner core material 

Beryllium copper 

Outer spherical shell material 

Aluminum 

Semimajor Axis 

12265 km 

Eccentricity 

0.004 

Inclination 

109.8° 

Perigee altitude 

5858 km 

Apogee altitude 

5958 km 

Perigee rate 

-0.214°/day 

Orbital period 

225.3 minutes 

Node rate 

+0.3437day 

Semimajor axis decay rate 

-1.1 mm/day 


LAGEOS is in a highly stable, nearly circular Earth orbit, at an altitude of 
approximately 5900 km. This high altitude, along with the satellite’s simple spherical 
shape and high mass density, minimize the effects of drag, solar radiation pressure, and 
other nonconservative forces. The effects of the short-wavelength geopotential 
coefficients, which are not known to the accuracy of the long- wavelength coefficients, 
are also minimized in this high altitude orbit. Thus, the motion of LAGEOS is quite 
predictable based on current satellite dynamical models. With both highly accurate 
range observations and a very predictable trajectory, LAGEOS is ideal for determining 
tracking station positions and distances between stations, geopotential coefficients, tidal 
amplitudes, nonconservative force parameters, polar motion parameters. Earth rotation 
parameters, tectonic plate motions, and other model parameters. Moreover, the 
temporal variations of such solutions can be determined more accurately and precisely 
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than from other satellites not specifically designed for precise orbit determination 
purposes. Many of the geodetic solutions obtained from LAGEOS data have also been 
verified by independent very long baseline interferometry (VLBI) and GPS techniques 
[Kolenkiewicz et al., 1985; Himwich et al., 1993; Dickey, 1993; McCarthy, 1993]. 

3.4 LAGEOS Geodetic Solutions 

Advances in SLR techniques have led to improvements in dynamic satellite 
geodesy methods and geodynamic solutions. By observing a satellite using SLR 
methods, its motion as it moves through the Earth’s geopotential field can be 
monitored. Thus, by modeling the satellite’s orbital motion due to gravitational and 
non-gravitational forces, SLR observations provide a link or connection between the 
satellite and the tracking stations observing the satellite. The parameters describing the 
satellite’s motion, the locations of the tracking stations, and other parameters in the 
model are estimated by minimizing the difference between the model orbit and the actual 
orbit in a least squares fashion through the use of the SLR observations. 

LAGEOS, being designed specifically for satellite geodesy and geodynamics 
research, has provided improved estimates of tracking stations and their motions since 
its mission began in 1976. One of the most current solutions is reported by Tapley et 
al. [1993]. Improved estimates of other geophysical parameters and their time 
variations have been made using LAGEOS SLR data as well. Lerch et al. [1985], 
Marsh et al. [1988], Marsh et al. [1990], Nerern et al. [1994], and others have detailed 
improvements in the Earth’s gravitational field from LAGEOS data. Tidal parameters 
[Christodoulidis et al., 1986], Earth rotation and polar motion parameters [Tapley et 
al., 1985; Pavlis et al., 1988; Caporali et al., 1990; Tapley et al., 1993], and secular 
variations in J 2 [Yoder et al., 1983; Rubincam, 1984; Cheng et al., 1989; Gegout et al., 
1991] have also been estimated based on LAGEOS data. Rubincam [1984] and others 
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have used LAGEOS data to study postglacial rebound. The estimates made from 
LAGEOS SLR observations are more meaningful than most other data types due to the 
accuracy of both the observations and the predicted orbit, which translates to improved 
accuracy of the estimates. The utilization of LAGEOS for satellite geodesy and 
geodynamics research has led to many advancements and achievements which would 
not have occurred as quickly without the benefit of LAGEOS. 

3.4.1 LAGEOS J 2 and J 3 Solutions 

Temporal variations in the low degree zonal harmonics of the Earth’s 
gravitational field are observable using SLR [Shum et al., 1987; Cheng et al., 1989, 

1990] . Some examples of the geophysical processes contributing to these variations 
are luni-solar tides, nontidal variations in the distribution of atmospheric mass and 
ocean water mass, meteorological mass redistribution, and postglacial rebound of the 
solid Earth [Chao et al., 1987; Cheng et al., 1989; Chao and Au, 1991]. Continental 
water storage and snow cover/loading can cause variations in the gravitational field also 
[Chao et al., 1987; Chao and O’Connor, 1988]. The seasonal gravitational variations 
(annual and semiannual), however, are dominant and are mainly due to solar influences 
on the Earth’s mass distribution. The primary solar influences include the seasonal 
redistribution of atmospheric mass, the seasonal redistribution of hydrospheric mass, 
and solar annual and solar semiannual solid Earth and ocean tides [Chao and Au, 

1991] , 

An example of a current solution for the temporal variations of the J 2 and J 3 
gravitational coefficients based on LAGEOS SLR data is that of Nerem et al. [1993]. 
As they detail, monthly estimates of the nontidal, nonsecular J 2 and J 3 variations are 
computed for the period from 1980 to 1990. Figure 3.2 presents the monthly estimates 
for the J 2 variation and Figure 3.3 shows the monthly estimates for the J 3 variation that 
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they computed. These solutions compute a constant J 2 and J 3 variation for every month 
during the 10 year period. Using the terminology described in the previous chapter, 
these are referred to as monthly boxcars. A resolution of one month is typical for 
solutions of temporal variations in gravitational coefficients based on SLR data. 
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Figure 3.2. Monthly Variations in the J 2 Gravitational Coefficient Computed from 
LAGEOS SLR Data from Nerem etal [1993]. 
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Figure 3.3. Monthly Variations in the J 3 Gravitational Coefficient Computed from 
LAGEOS SLR Data from Nerem et al. [1993]. 


3.4.2 LAGEOS Along Track Drag Solutions 

After subtracting most of the known forces acting on LAGEOS, there is still an 
along track deceleration (drag) which has been observed and reduces the semimajor 
axis by about 1 mm per day (see Tapley et al [1993] for a more detailed look at the 
long term estimates of this drag and an interesting discussion on the drag in general). 
Previous studies have addressed the possible origins and proposed models for this drag 
[Rubincam, 1982, 1987, 1988, 1990], which are mostly thermal in nature. 

A current solution for the anomalous along track drag on LAGEOS is computed 
by Tapley etal [1993]. This along track drag remains after most of the known forces 
acting on LAGEOS are taken into account. Estimates of this along track drag are 
necessary when estimating other geodynamic or geodetic parameters from LAGEOS 
SLR data to prevent this acceleration from aliasing into estimates. Typically, the along 
track drag is estimated every 15 days (15 day boxcars). Figure 3.4 shows the boxcar 
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estimates computed by Tapley et al. based on LAGEOS data from May 30, 1976 to 
May 30, 1993. The along track drag has a mean of -3.5 picometer/s 2 over the arc. 
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Figure 3.4. 15 Day Estimates of the Along Track Acceleration for LAGEOS from 

Tapley etal [1993]. 


3.5 Proposed Benefits of Process Noise Filtering for LAGEOS 

Most LAGEOS data analysis techniques determine geodetic and geodynamic 
parameters by solving for them as if they were constants over the specific time interval 
of the data arc through the use of a standard filter. For example. Smith et al [1985], 
estimate polar position values and variations in universal time every 5 days, and 
tracking station locations and other parameters every 30 days. Smith et al [1991] 
estimate parameters as constants over fixed time intervals as well. Aside from being 
linked from one time interval to the next by the satellite orbital parameters, these 
constant parameters are effectively decoupled from one time interval to the next. In 
other words, the estimates for the constants over the time interval are primarily based 
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on the observations made during the time interval Estimates of purely kinematic 
parameters are based entirely on observations made during the particular time interval, 
while estimates of dynamic parameters are still linked from one time interval to the next 
by the satellite orbit. If these parameters, in reality, do not change over the given time 
interval, then model error is not introduced into the solution. If, however, these 
parameters do vary over the given time interval, then model error will corrupt the 
solution. By limiting the time interval such that the given parameters do not change in 
any measurable way over that interval, then model error can be kept at a minimum. A 
time series of a specific parameter can then be constructed by joining the consecutive 
fixed interval estimates of the parameter (boxcars). The fundamental drawback to 
generating temporal solutions in this manner is the actuality that estimates will degrade 
if the time interval for the boxcar estimate is either too short or too long. If the time 
interval is too short, then accurate estimates may not be possible given the available data 
or the observability of the variation. If the time interval is too long, then the parameter 
may vary significantly from the constant estimate over the interval, thus resulting in 
adverse effects from model error. The shorter boxcar time intervals may be desirable 
from a resolution perspective, while the longer boxcar time intervals may help to reduce 
the effects of measurement noise and erroneous data. 

By estimating these parameters stochastically using a process noise filter, long 
data arcs may be used to generate a solution from LAGEOS SLR data without the 
division of the arc into smaller arcs for the estimates of known time varying parameters. 
The strength of the solution may be enhanced by using a long time interval (the length 
of the entire arc for each parameter), and the temporal resolution of the solution of the 
parameters may also improve. Along with the entire geodynamic solution, the estimate 
of the satellite orbit over time may be improved as well. The time variations of 
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parameters are found by using the entire long arc of data, even though the parameters 
may change significantly over the long interval. For example, variations in gravity field 
coefficients or variations in tracking station positions could be estimated with improved 
resolution by using long arcs spanning years to decades rather than the same long arcs 
divided into shorter monthly time intervals for these parameters. An estimate for each 
process noise parameter is computed at every time that an observation is made over the 
data arc. Thus, the estimate at each individual time for a process noise parameter is 
based on the observations made over a rough, adjustable time interval centered on the 
time of the estimate. The length of this rough time interval, or correlation window, is 
determined by a time constant T, for the particular parameter. Effectively, observations 
made outside of this correlation window in time have little or no effect on the estimate 
made at the time at the center of the window. Parameters which are highly correlated in 
time (slowly varying) may have longer correlation windows than parameters which 
vary more rapidly in time. Depending on the observability of the particular variation in 
the observation, a parameter may be estimated more accurately as a process noise 
parameter. For variations which are not effectively observed in the satellite 
observations, the process noise estimates may be no more accurate than standard, 
consecutive fixed interval estimates. The possible benefits of estimating particular 
parameters stochastically based on LAGEOS SLR data are assessed in this study. 



CHAPTER 4 


SIMULATION MODEL 


4.1 Introduction 

Model errors of some type exist within any real satellite observation. This is a 
result of the fact that the actual forces causing the motion of the satellite differ from the 
forces that can be represented in a mathematical force model. One can not know the 
true variations from the dynamical force model. This results in an imperfect knowledge 
of the true trajectory of a satellite if real observations are used. The best one can do is 
determine the best estimate of the trajectory, which is limited by how well your model 
represents the true trajectory. 

To fully understand how well model errors can be resolved by estimating 
parameters stochastically, the model errors themselves must be known. In this study, 
the true temporal variations of the particular parameters being estimated must be known 
perfectly before conclusions can be drawn as to the effectiveness of particular filtering 
methods in estimating these variations. Thus, simulated LAGEOS SLR observations 
are generated with specific temporal variations of parameters (referred to hereafter as 
model deviation signals or model signals) built into the model used in generating the 
observations. This permits the direct comparison of the estimates to the known truth. 

LAGEOS SLR data is generated for a period of one year for this study. 
Solutions obtained by processing this simulated one year arc both with stochastic and 
non-stochastic filtering techniques are compared. The process by which this data is 
generated is now detailed. A three year LAGEOS SLR data arc, processed only by the 
stochastic filter, is also generated following the same simulation procedure. The laser 
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ranging measurements are computed based on a specified tracking network of laser 
ranging stations and a satellite in the same orbit as LAGEOS, acting under a specified 
dynamical model. This chapter details the models used and the process by which the 
simulated observations are generated. 

4.2 Measurement Model 

The measurement used in the simulation is a range observable, which 
corresponds to a typical SLR measurement. The magnitude of the difference between 
the satellite position vector and the tracking station position vector is determined: 

Z = llr s -r, II (4.1) 

where Z is the SLR measurement (range), r s is the satellite position vector, and r t is the 
tracking station position vector. This measurement simply represents the distance 
(range) from the tracking station to the satellite. Figure 4.1 shows the geometry of the 
SLR measurement. The actual determination of real laser range measurements is a bit 
more complex, as it involves corrections for atmospheric refraction, instrument delays, 
relativistic effects, and the finite speed of light. 
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The final, corrected laser range measurement is what is computed geometrically 
in equation 4.1 and used in this simulation. Finally, white noise with 1 cm RMS is 
added to all measurements once generated. This random error corresponds to the 
current ideal levels of accuracy and precision associated with SLR observations. Actual 
data may suffer from biases which are being ignored in this simulation. 

4.2.1 Tracking Network 

Typically, dozens of laser ranging tracking stations are able to track LAGEOS 
using the SLR technique. Over 100 tracking station sites around the world have made 
laser ranges to LAGEOS since its launch. In generating the simulated observations for 
this study, only eight tracking stations are used. This conservative tracking network 
was chosen to represent a worst case tracking scenario. Also, using a subset of the 
actual tracking network will reveal any problems that might result from the lack of a 
dense data distribution. 
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Figure 4.2 shows the tracking station network. Six of the eight tracking 
stations are located in the northern hemisphere. This reflects the fact that the majority 
of SLR tracking stations are in the northern hemisphere. 



Figure 4.2. Tracking Station Network for LAGEOS Data Simulation With 20 p 
Visibility Masks Shown. 

The tracking stations operate from 6:00 PM to 6:00 AM local time. Four of the 
northern hemisphere stations track five days per week (Monday through Friday), as is 
the case with many NASA stations, while the rest of the stations track every day. 

Observations are generated every three minutes (representing compressed 
normal point observations) for a specific tracking station only if the tracking station is 
operating during the particular pass, and the satellite is above 20° elevation. All 
observations are decimated by randomly eliminating 75% of the passes in an attempt to 
model data outages due to weather problems or inoperative tracking stations. Figure 
4.3a shows the histogram of the simulated LAGEOS data over the one year simulation 
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period and Figure 4.3b shows a histogram of actual LAGEOS data over the same time 
period. 
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Figure 4.3a. Histogram of Simulated LAGEOS Laser Ranging Data. 
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Figure 4.3b. Histogram of Actual LAGEOS Laser Ranging Data. 
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The simulation generates a conservative amount of data relative to the actual 
amount of LAGEOS data that is available from the same time period. The number of 
passes in each six days during the simulated LAGEOS arc is approximately one-half the 
number of passes in each six days during the actual LAGEOS arc. This ensures that 
any conclusions drawn based on the results of this simulation will not be improperly 
due to excessive or unrealistic data density. Table 4.1 shows the number of passes and 
observations for each of the tracking stations during the one year simulation. The days 
of the week that the stations operate is also summarized. 


Table 4.1. Tracking Stations and Number of Passes and Observations Generated 
During the One Year Simulation. Station operation schedule is also summarized. 


Tracking 

Station 

Station 

ID 

Number of 
Passes 

Number of 
Observations 

Station Operation 
Schedule 

GRF105 
(GSFC, USA) 

7105 

151 

1897 

Monday-Friday 

QUINC2 
(Quincy, USA) 

7109 

143 

1777 

Monday-Friday 

HOLLAS 
(Maui, USA) 

7210 

118 

1403 

Monday-Friday 

WETZEL 

(Wettzell, Germany) 

7834 

135 

1765 

Monday-Friday 

ARELAS 
(Arequipa, USA) 

7907 

186 

2261 

Monday-Sunday 

MATERA 
(Matera, Italy) 

7939 

208 

2712 

Monday-Sunday 

SHO 

(Simosato, Japan) 

7838 

201 

2491 

Monday-Sunday 

YARAG 

(Yaragadee, Australia) 

7090 

238 

2586 

Monday-Sunday 

Totals 


1380 

16892 
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4.3 Dynamical Model 

A simplified dynamical force model is used in this study. The gravitational 
forces are modeled by a central body term (fi, the Earth’s gravitational coefficient) and 
the zonal nonspherical geopotential coefficients up to degree five (7 2 , 7 3 , 7 4 , and J s ). 
The reference value for fi is taken from Ries et al. [1992], and the zonal coefficients are 
taken from the JGM-2 [Nerem et al., 1994] gravity model. Thus, the geopotential 
model is longitudinally symmetric. Table 4.2 summarizes the geopotential model used. 

Table 4.2. Earth’s Gravitational Coefficient and JGM-2 Coefficients (Unnormalized) 
used in Geopotential Model. 


Coefficient 

Nominal Value 

Standard Deviation (o) 

V 

398600.4415 kmVs 2 

0.0008 km 3 /s 2 

h 

1082627.0 x 10' 9 

0.244 x 10" 9 

7j 

-2532.308 x 10' 9 

0.690 x 10~ 9 

7 4 

-1620.430 x 10' 9 

0.780 x 10~ 9 

7 5 

-227.0711 x 10“ 9 

0.521 x 10~ 9 


The anomalous along track drag observed in the LAGEOS orbit is an ideal 
parameter to estimate stochastically since the forces causing the drag are not completely 
understood. In this study, this particular along track drag is considered the total drag. 
An empirical drag model is used to model this acceleration. The empirical function 
used is 


r 

Ilf II 


D=C. 


(4.2) 
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where C, is the along track parameter for the drag acceleration. The nominal average 
value for C t used in this study is the mean of the observed value, or -3.5 x 10“ 12 m/s 2 
(-3.5 picometer/s 2 ) [Tapley et al., 1993]. 

The satellite state deviation vector, x, that is used in the simulation is 

x = [ C t J 2 J 3 Js, Js /l x y z x y z ] T 

It is noted that estimates of this state vector are state deviations from the reference state 
based on the dynamical force model. 

4.4 True Model Deviation Signals Added to Model Parameters 

The dynamical model parameters described above remain fixed in the model 
throughout the filtering process. However, specific temporal variations, or model 
deviation signals, are introduced into the C„ J 2 , J 3 , J 4 , and J 5 model parameters in the 
dynamical force model as the simulated SLR measurements are generated. The 
temporal estimates of the these signals after filtering may then be compared to the 
known true model deviation signals present in the data. Again, these model deviation 
signals are not modeled in the dynamical model used in the filtering of the data. They 
are only added to the constant nominal values of the respective model parameter during 
the generation of the simulated data. 

Some basic assumptions are made in defining these realistic model deviation 
signals. In general, the temporal signals are based on previously reported estimates or 
models of the temporal variations for the particular parameters. The true model 
deviation signals for C t , J 2 , and / 3 are based on previously reported estimates. While 
unknown, it is presumed that these parameters vary in some continuous manner. The 
previously reported estimates for these signals are then interpolated using a natural 
cubic spline in order to produce a smooth, continuous signal to use for the truth in this 
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study. The natural cubic spline interpolation forces the interpolation through the 
support points (previously reported estimates) while generating a smooth function. The 
true model deviation signals for J 4 and J 5 are based on proposed models. The true 
signal used for these parameters is inherently continuous since it is based on a model. 
No assumptions are made as to the expected averages (biases) of the model deviation 
signals for each of the parameters. While there is no expectation that the averages for 
these model deviation signals will be zero, no separate estimate is made for these 
biases. That is, the total model deviation signal for each parameter is estimated as a 
single deviation in the filtering process rather than a bias plus a variation from the bias. 
This approach should be valid as long as the biases of the model deviation signals 
remain on the order of magnitude of the variations from the bias, which is the 
circumstance for the model deviation signals used in this study and detailed below. 

All of the model deviation signals for the one year arc are based on estimates 
made during the one year period from January 1, 1986 to December 31, 1986. The 
true model deviation signals introduced into the simulated data are detailed below. 

4.4.1 Drag Model Deviation Signal 

For the along track drag parameter C t , the model deviation signal shown in 
Figure 4.4 is added to the nominal C t value of -3.5 picometer/s 2 during the simulation 


of the observations. 
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Figure 4.4. True Model Deviation Signal for the Q Parameter Based on Current Best 
Estimates. 


This model deviation signal is taken from Tapley et al [1993] which gives 15 day 
estimates for the observed along track drag for LAGEOS over a 14 year period. The 15 
day estimates falling between January 1, 1986 and December 31, 1986 are also shown 
in Figure 4.4 for reference. This one year of interpolated along track drag variation, 
with the average value subtracted out, is what is shown in Figure 4.4. 


4.4.2 J 2 and J 3 Model Deviation Signal 

The model deviation signal for the J 2 coefficient consists of a secular and 
nonsecular term. The model deviation signal for the J 3 coefficientis purely nonsecular. 
The secular rate used for J 2 ( j 2 ) is -2.6 x 10' n /yr [Nerem et al, 1993]. The total 
model deviation signal for each coefficientrepresents the nontidal temporal variations in 
J 2 and J 3 . 

The nonsecular part of the J 2 and J 3 model deviation signal is taken from Nerem 
et al. [1993] which gives monthly estimates for the nonsecular variations of J 2 and J 3 
over the time period from 1980 to 1989. The estimates used are those they computed 
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from atmospheric pressure data (with no correction for the inverted barometer effect). 
These estimates are used in this study since they represent the current best estimates of 
the true J 2 and J 3 nonsecular variations. (Similar estimates of the J 2 and J 3 nonsecular 
variations exist based on LAGEOS SLR data, but are “effective” estimates as they 
include nonsecular variations from higher degree zonals which cannot be separated 
without independent data from satellites in different orbits.) Their monthly estimates 
for the nonsecular variations from January, 1986 to December, 1986 are interpolated 
with a natural cubic spline. This one year of interpolated J 2 and J 3 nonsecular variation, 
with the secular variation for J 2 added back, is what is shown in Figures 4.5 and 4.6 
with the monthly estimates they are based on shown for reference (unnormalized). 
This is the total model deviation signal (total temporal variation) added to the nominal 
JGM-2 J 2 and / 3 values. 



Figure 4.5. True Model Deviation Signal for the J 2 Parameter Based on Current Best 
Estimates. 
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Figure 4.6. True Model Deviation Signal for the J 3 Parameter Based on Current Best 
Estimates. 


4.4.3 ,/ 4 and J s Model Deviation Signal 

For the / 4 and i 5 coefficients, the model deviation signal shown in Figures 4.7 
and 4.8 is added to the nominal JGM-2 7 4 and J 5 values (unnormalized). It should be 
noted that these model deviation signals are only included in part of this study, as 
discussed in chapter 5. 
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Figure 4.7. True Model Deviation Signal for the J 4 Parameter Based on Current 
Models. 



Figure 4.8. True Model Deviation Signal for the J 5 Parameter Based on Current 
Models. 


These J 4 and J 5 model deviation signals are taken from Chao and Au [1991] 
which gives the amplitude and phase of the seasonal variations from the mean of J 4 and 
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J 5 (among other coefficients) based on global surface pressure data from 1980 to 1988. 
Table 4.3 details the characteristics which are pertinent to this study. 


Table 4.3. Amplitude and Phase of Seasonal Variations in J 4 and J 5 due to Atmospheric 
Mass Redistribution Without the Oceanic Inverted Barometer Effect (from Chao and Au 
[1991]). Phase with respect to the sine convention with t = 0 on January 1. 



Annual Variation J 

Semiannual Variation 

Coefficient 

Amplitude 

Phase 

Amplitude 

Phase 

h 

1.52 x 10" 10 

-155° 

1.12 x 1(T 10 

-101° 

h 

1.38 x 1(T 10 

o 

T-H 

1 

0.44 x 10 10 

62° 


For this study, the temporal variations from January 1, 1986 to December 31, 1986 are 
computed from a function based on the amplitudes and phases noted above. Also, 
biases of 7.80 x 10“ 10 and 5.21 x 10~ 10 are included in the model deviation signals for J 4 
and J 5 respectively. These biases represent the JGM-2 standard deviations (o) for the 
reference values of 7 4 and J 5 that are used in the dynamical model. This total model 
deviation signal (bias and temporal variation) is what is shown in Figures 4.7 and 4.8. 

4.5 Summary of Data Simulation 

To summarize, the satellite trajectory is numerically integrated using the 
dynamical model and the model deviation signals. These model deviation signals are 
continuous and based on a natural cubic spline interpolation of 15 day current estimates 
for C, and one month estimates J 2 and J 3 for the period from 1 January 1986 to 1 
December 1986. The continuous model deviation signals for J 4 and J s are based on 
current models. Range data from specified tracking stations to this trajectory is then 
computed at three minute intervals (based on observability and station operation criteria) 
and saved with 1 cm RMS random noise added to each range measurement. This data, 
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comprised of a time from epoch, tracking station number, and range measurement, 
makes up the simulated SLR data set used in this analysis. 



CHAPTER 5 


FILTERING RESULTS 


5.1 Introduction 

This chapter summarizes the solutions obtained from both the standard and 
stochastic orbit determination filtering modes. Both methods are used independently to 
determine solutions from the same set of simulated observations. A simulated one year 
data arc is processed both with and without the 7 4 and 7 5 model deviation signals 
present in the data. Also, a simulated three year data arc (without the 7, or J s model 
deviation signals present) is processed using the stochastic filtering mode and the 
resulting solution is presented. 

5.2 One Year Arc: C t , J 2 , and J 3 Model Deviation Signals 

A simulated one year data arc with model deviation signals present in C„ J 2 , and 
7 3 and no model deviation signals present in 7 4 or 7 5 is processed. For both the 
standard and stochastic filtering methods, the one year data arc is processed with a 
single estimate made for the satellite state ( r and r ), and no correction made for fi, 7 4 , 
or J 5 (since no model deviation signal is present in the fi, 7 4 , or 7 5 parameters, no 
correction is needed). In section 5.3, data that contains 7 4 and 7 5 model deviation 
signals is processed. For the standard SRIF mode, two iterations are performed in 
generating the solution. The reference trajectory for the second iteration uses the 
corrections from the first iteration. The estimates for C„ J 2 , and 7 3 from the first 
iteration are used in the dynamical model during the second iteration as well. For the 
process noise SRIF mode, a preliminary correction is made only to the satellite state 
while fixing all other parameters. This gives a reference orbit that is sufficiently close 
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to the true orbit over the entire arc, thus minimizing linearization errors. A second 
iteration is performed using the process noise mode once this preliminary, non- 
stochastic correction is made. 

5.2.1 Standard SRIF Using Boxcars 

The standard SRIF is implemented using consecutive 15 day estimates for C t 
throughout the arc. That is, a single, constant C t estimate is made for the first 15 days 
of the arc, a second, constant C t estimate for the second 15 days of the arc, and so on. 
This results in 24 consecutive C t estimates over the one year arc. Thus, the solution for 
the C t variation over the one year period is defined by joining these 24 C t “boxcars” 
together. Again, the term boxcar is used henceforth to denote this type of solution 
based on its appearance as the discontinuous, constant estimates are consecutively 
joined together. Similarly, one month estimates are made for J 2 and J 3 throughout the 
one year arc. This results in 12 consecutive boxcars for the J 2 and / 3 solution over the 
one year arc. Thus, the simultaneous solution for the 24 C, boxcars, the 12 J 2 boxcars, 
the 12 7 3 boxcars, and the satellite position and velocity is required. 

First, the boxcar estimates for the C„ J 2 , and / 3 temporal variations are 
compared to the truth. Figure 5.1 shows the true C, and the estimated C t and Figure 
5.2 shows the C, residuals, or the true C t minus the estimated C t . 
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Figure 5.1. True C t and Estimated Q Using Boxcars. 



Figure 5.2. C, Residuals Using Boxcars. 


The boxcar estimates for Q track the true signal quite well. This is consistent with the 
success that others have had using this approach to estimate C, variations with actual 
LAGEOS data [Tapley et aL, 1993]. Figures 5.3 and 5.4 show the true J 2 and 
estimated J 2 and the J 2 residuals respectively. 
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Months Past 1 January 1986 
Figure 5.3. True J 2 and Estimated J 2 Using Boxcars. 



Months Past 1 January 1986 


Figure 5.4. J 2 Residuals Using Boxcars. 


Likewise, Figures 5.5 and 5.6 show the true J 3 and estimated J 3 and the J 3 residuals 
respectively. 
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Figure 5.5. True J 3 and Estimated J 3 Using Boxcars. 



Figure 5.6. J 3 Residuals Using Boxcars. 

Again, the boxcar estimates track the true signals quite well for both the J 2 and J 3 
parameters. This too is consistent with the success that others have had using this 
approach to estimate J 2 and J 3 variations with actual LAGEOS data [Nerem et aL, 
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1993]. Table 5.1 summarizes the RMS statistics for the C„ J 2 , and / 3 residuals. The 
RMS statistics of the respective model deviation signals are also shown for reference. 

Table 5.1. RMS for C„ J 2 , and J 3 Residuals Using Boxcars. 


Parameter 

Model Deviation 
Signal RMS 

Residual RMS 
Using Boxcars 

c t 

1.37 x 1(T 12 m/s 2 

2.21 x 10 13 m/s 2 

J 2 

2.85 x KT 10 

8.66 x 10-' 1 

A 

3.10 x KT 10 

9.37 x lO - " 


Next, the differences between the true orbit and the estimated orbit are analyzed 
from a positional standpoint. The radial, transverse, and normal (RTN) residuals are 
shown in Figures 5.7 to 5.9. These residuals represent the true position minus the 
estimated position in each direction respectively. The radial direction is the direction of 
the satellite position vector r. The normal direction is the direction resulting from the 
matrix cross product of the satellite position vector r and the satellite velocity vector r 
(the normal direction is perpendicular to the satellite’s plane of motion). The transverse 
direction results from the matrix cross product of the normal and radial directions. 






Transverse Residuals (cm) « Radial Residuals (cm) 


J] 



Figure 5.8. Transverse Residuals Using Boxcars. 





5.10. This is simply the magnitude of the positional difference between the hue 
position and the estimated position of the satellite at each time. 



Figure 5.10. 
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Figures A. 1 through A.6 in Appendix A show the orbit residuals in terms of the 
Keplerian elements (semi-major axis a , eccentricity e, inclination /, argument of 

periapse co, longitude of ascending node Q, and argument of latitude co +/, where / is 
true anomaly). 

Finally, the range residuals (observed range minus computed range) are 
determined. Figure 5.11 shows the range residuals, and Table 5.2 summarizes the 
RMS statistics for the RTN, 3-d position, and range residuals. Since 1 cm RMS noise 
is present in the simulated observations, it is expected that the RMS of the range 
residuals will approach 1 cm as the estimated solution approaches the truth. 



0 1 23456789 10 11 12 

Months Past 1 January 1986 


Figure 5.11. Range Residuals Using Boxcars. 
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Table 5.2. RMS for RTN and Range Residuals Using Boxcars. 


Residual 

RMS (cm) 

Radial 

0.68 

Transverse 

6.43 

Normal 

3.61 

3-d Position 

7.40 

Range 

2.26 


These residuals, along with the individual C t , J 2 , and J 3 parameter residuals, form the 
basis for the statistical comparison of the standard SRIF boxcar solution to the process 
noise SRIF solution. 

5.2.2 Process Noise SRIF 

The process noise SRIF is implemented with stochastic estimates made for C t , 
J 2 , and J 3 . Stochastic estimates of each parameter are made at every time that an 
observation exists. Specific values for rand <7 (as defined in chapter 2) are chosen for 
each parameter based on the expected time correlation and amplitude of the model 
deviation signal for that parameter. Table 5.3 summarizes the values used for each 
parameter. 
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Table 5.3. Values for the Time Correlation Constant t, and Steady State Standard 
Deviation a, for One Year Arc with C„ J 2 , and J 3 Model Deviation Signals. 


Parameter 

t (years) 

a 

c, 

1 

12 

3.5 x 1(T 12 m/s 2 

J 2 

1 

4 

1.5 x 10 -9 

h 

1 

4 

3.0 x IQ' 9 


The stochastic estimates for the C„ J 2 , and J 3 temporal variations are compared 
to the truth. Figure 5.12 shows the true C, and the estimated C, and Figure 5.13 shows 
the C t residuals, or the true C t minus the estimated C t . 



Figure 5. 12. True C, and Estimated C x Using Process Noise. 
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Figure 5.13. C, Residuals Using Process Noise. 

The stochastic estimate of C t appears excellent. Relative to the previous boxcar estimate 
for Q, the stochastic estimate appears to have much better temporal resolution and 
accuracy. Figures 5.14 and 5.15 show the true J 2 and estimated J 2 and the J 2 residuals 
respectively. 



Figure 5.14. True J 2 and Estimated J 2 Using Process Noise. 
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Figure 5.15. J 2 Residuals Using Process Noise. 


Likewise, Figures 5.16 and 5.17 show the true/ 3 and estimated/} and the / 3 residuals 
respectively. 



Figure 5.16. True / 3 and Estimated/ 3 Using Process Noise. 
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Figure 5.17. i 3 Residuals Using Process Noise. 

The stochastic estimates of J 2 and J 3 are also excellent. Relative to the previous boxcar 
estimates for J 2 and 7 3 , the stochastic estimates appear to have much better temporal 
resolution and accuracy. Table 5.4 summarizes the RMS statistics for the C„ J 2 , and / 3 
residuals with the respective RMS statistics from the previous boxcar solution shown 
for comparison (from Table 5.1). The RMS statistics of the respective model deviation 
signals are also shown for reference. 


Table 5.4. RMS for C t , J 2 , and J 3 Residuals Using Process Noise and Boxcars. 


Parameter 

Model Deviation 
Signal RMS 

Residual RMS 
Using Boxcars 

Residual RMS 
Using Process Noise 

Q 

1.37 x 10" 12 m/s 2 

2.21 x 10 -13 m/s 2 

5.91 x 10~ 14 m/s 2 

J 2 

2.85 x 1(T 10 

8.66 x 10 _u 

1.96 x 10-" 

J 3 

3.10 x lO' 10 

9.37 x 10“" 

2.36 x 10-" 
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The process noise residuals for C„ J 2 , and J 3 are noticeably improved relative to the 
respective boxcar residuals. 

Next, the RTN positional residuals are shown in Figures 5.18 to 5.20. The 
total RTN, or 3-d position residuals are shown in Figure 5.21. 



0 1 23456789 10 11 12 


Months Past 1 January 1986 
Figure 5.18. Radial Residuals Using Process Noise. 



0 1 23456789 10 11 12 


Months Past 1 January 1986 

Figure 5.19. Transverse Residuals Using Process Noise. 





3-d Position Residuals (cm) 01 Normal Residuals (cm) 
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Months Past 1 January 1986 
.20. Normal Residuals Using Process Noise. 
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0 1 23456789 10 11 12 

Months Past 1 January 1986 

Figure 5.21. 3-d Position Residuals Using Process Noise. 




Figures B.l through B.6 in Appendix B show the orbit residuals in terms of the 
Keplerian elements (semi-major axis a, eccentricity e, inclination i, argument of 
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periapse co, longitude of ascending node Q, and argument of latitude co +/, where / is 
true anomaly). 

Figure 5.22 shows the range residuals, and Table 5.5 summarizes the RMS 
statistics for the RTN, 3-d position, and range residuals (with comparison to the 
respective boxcar statistics from Table 5.2 for reference). Again, since 1 cm RMS 
noise is present in the simulated observations, it is expected that the RMS of the range 
residuals will approach 1 cm as the estimated solution approaches the truth. 



Months Past 1 January 1986 


Figure 5.22. Range Residuals Using Process Noise. 
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Table 5.5. RMS for RTN and Range Residuals Using Process Noise and Boxcars. 


Residual 

RMS Using 
Boxcars (cm) 

RMS Using 
Process Noise (cm) 

Radial 

0.68 

0.47 

Transverse 

6.43 

0.52 

Normal 

3.61 

0.80 

3-d Position 

7.40 

1.07 

Range 

2.26 

1.11 


Figures 5.23 to 5.27 show the RTN, 3-d position, and range residuals from 
both the standard boxcar SRIF and the process noise SRIF solutions side by side on 
the same scale for visual comparison. 




(a) (b) 

Figure 5.23. Radial Residuals Using Boxcars (a) and Process Noise (b). 
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(a) (b) 

Figure 5.24. Transverse Residuals Using Boxcars (a) and Process Noise (b). 



(a) 



0 


12 


Months Past 1 January 1986 
(b) 


Figure 5.25. Normal Residuals Using Boxcars (a) and Process Noise (b). 
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(a) (b) 

Figure 5.26. 3-d Position Residuals Using Boxcars (a) and Process Noise (b). 



0369 12 0369 12 

Months Past 1 January 1986 Months Past 1 January 1986 


(a) (b) 

Figure 5.27. Range Residuals Using Boxcars (a) and Process Noise (b). 

Clearly, the stochastic filtering method produces a temporal solution for Q, J 2 , 
and / 3 which is more accurate than the standard boxcar solution. The RMS for the C t , 
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J 2 , and J 3 residuals are reduced to approximately 25% of the respective residuals using 
the standard boxcar method (a 75% improvement). Similarly, by using the stochastic 
filtering method the total 3-d positional residual is reduced to approximately 15% of the 
3-d positional residual resulting from the standard boxcar filtering method (an 85% 
improvement). The range residuals are reduced by about 50%. 

5.3 One Year Arc: C t , J 2 ,J 3 , J 4 , and J 5 Model Deviation Signals 

Next, a simulated one year data arc with model deviation signals present for C„ 
J 2 , /?, J A , and J 5 is processed. Again, the one year data arc is processed with a single 
estimate made for the satellite state, and no estimate made for /i, J A , or 7 5 , even though 
J 4 and J s model deviation signals are present in the data for this case. This 
configuration assesses the ability of each filter to resolve the C, variation and the 
effective J 2 and J 3 temporal variations since J 4 and J$ model deviation signals are present 
but not estimated (a discussion on the effective J 2 and J 3 temporal variations will follow 
shortly). This situation more closely parallels typical filtering scenarios which estimate 
temporal variations in particular parameters in the presence of other parameter 
variations/model deviation signals which are not estimated. Aside from the introduction 
of the 7 4 and J 5 model deviation signal, the analysis parallels the previous scenario. 

5.3.1 Standard SRIF Using Boxcars 

The standard SRIF is implemented using 24 consecutive 15 day estimates for C, 
and 12 consecutive one month estimates for J 2 and J 3 throughout the arc. First, the 
boxcar estimates for the C„ J 2 , and J 3 temporal variations are compared to the truth. 
Figure 5.28 shows the true C t and the estimated C, and Figure 5.29 shows the C, 
residuals, or the true C, minus the estimated C t . The boxcar estimates for C, appear to 
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be as accurate as in the previous boxcar case where no J 4 or J 5 model deviation signals 
are present. 



Figure 5.28. True C, and Estimated Q Using Boxcars with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 



Figure 5.29. C t Residuals Using Boxcars with J 4 and J 5 Model Deviation Signals 
Present but not Estimated. 
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Figures 5.30 and 5.31 show the true and estimated J 2 and true and estimated J 3 
respectively. 



Figure 5.30. True7 2 and Estimated / 2 Using Boxcars with J A and J 5 Model Deviation 
Signals Present but not Estimated. 



Figure 5.31. True J 3 and Estimated J 3 Using Boxcars with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 




79 


It is clear from Figures 5.30 and 5.31 that the estimates of the J 2 and J 3 
variations with J 4 and J s model deviation signals present but not estimated do not 
correspond to the true J 2 and i 3 model deviation signals very well. A bias differentiates 
the two. This is a result of the fact that the estimates of the J 2 and 7 3 variations are 
estimates of the “effective” J 2 and / 3 variations. A discussion of effective J 2 and J 3 
signals now follows. 

5.3.1.1 Effective J 2 and J 3 Model Deviation Signals 

For this part of the analysis, J 4 and J s model deviation signals are present in the 
data, but not estimated. In this scenario, estimates of the J 2 and / 3 variations are 
estimates of the “effective” J 2 and 7 3 variations. That is, the estimate of the J 2 variation 
is a value that alone would have the same effect (acceleration) on the satellite that the J 2 
and J 4 variations cause. Put another way, the J 2 and J 4 variations in the data are 
combined into an single (effective) J 2 variation. This effective J 2 variation is estimated 
by the filter since no estimate is being made for the J 4 variations. Further, if the J 4 
variations were estimated, the filter could not separate the J 2 and J 4 variations based on 
observations only from LAGEOS. The two variations are not separable without 
observations from other satellites. Likewise, the same relationship exists between J 3 
and / 5 . For the most part, when temporal variations of low degree zonal geopotential 
coefficients are estimated from single satellite data, they are effective estimates since 
variations in the higher degree zonals are not estimated separately. Thus, the effect of 
variations in all higher even degree zonals can be aliased into the variations of even low 
degree zonals that are estimated. Similarly, the effect of variations in all higher odd 
degree zonals can alias into the variations of odd low degree zonals that are estimated. 

The “true” effective J 2 (and true effective J 2 variation) is determined as follows: 
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The secular variation in the longitude of node Q, due to J 2 may be written 


dQ_ _ 
dt 



(5.1) 


where n is the satellite mean motion (nV = (i), a e is the semi-major axis of the central 
body’s reference ellipsoid, i is the satellite inclination, and a is the satellite semi-major 
axis. 

Likewise, the secular variation in Q. due to J A may be written 


= - — « - - - QS -- (l + ~r e 2 ) (— sin 2 i - l) 7 4 
dt 4 i- e 2 4 ' 2 /V4 > 


(5.2) 


The expression for the secular variation in Q, due to both J 2 and J 4 may be equated with 
an expression for the secular variation in Q due to an effective J 2 : 



(5.3) 


where is the effective J 2 . 


Thus, the effective J 2 (and effective J 2 variation) is simply a linear function of J 2 and J 4 : 



(5.4) 


For LAGEOS, the expression for the secular variation in the longitude of node due to J 2 
and J 4 (equations 5.1 and 5.2) is suitable for determining an effective J 2 since the 
LAGEOS orbit is very sensitive to secular changes in the longitude of node. For other 
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satellites, using expressions for the secular variation in argument of periapse, co , or 
mean anomaly M, may be preferable if they are not suitably sensitive to the node rate. 

The “true” effective 7 3 (and true effective 7 3 variation) is determined as follows: 
The long period variation in the eccentricity e, due to 7 3 may be written 



(5.5) 


where a is the satellite semi-major axis, a e is the semi-major axis of the central body’s 
reference ellipsoid, i is the satellite inclination, n is the satellite mean motion, and (O is 
the satellite argument of periapse. 

Likewise, the long period variation in eccentricity due to 7 5 may be written 



© 


5 / 4 + 3e 2 

l(l - 2 ) 4 


15 

16 


sin i + 


105 

32 


sin 3 i - sin 5 ij cos coJ 5 


(5.6) 


The expression for the long period variation in eccentricity due to both J 3 and J 5 may be 
equated with an expression for the long period variation in eccentricity due to an 
effective 7 3 : 



where is the effective J 3 . 


Q'Tu 
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Thus, the effective J 3 (and effective J 3 variation) for a particular satellite is simply a 
linear function of J 3 and J s : 


7L — + 




Again, for LAGEOS, the expression for the long period variation in eccentricity due to 
7 3 and J 5 (equations 5.5 and 5.6) is suitable for determining an effective J 3 since the 
LAGEOS orbit is sensitive to long period changes in eccentricity. For other satellites, 
using expressions for the long period variation in a different element may be preferable 
if they are not suitably sensitive to variations in eccentricity. 

Figure 5.32 shows the effective J 2 model deviation signal based on the actual J 2 
and J A model deviation signals shown in Figures 4.5 and 4.7. It is shown with the 
actual J 2 model deviation signal for reference. Thus, the actual J 2 model deviation 
signal combined with the actual J A model deviation signal produces the same 
acceleration to LAGEOS as the effective J 2 model deviation signal combined with no J A 
model deviation signal (just a constant reference J A ) . 
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Figure 5.32. Effective and Actual J 2 Model Deviation Signals. 

Figure 5.33 shows the effective J 3 model deviation signal based on the actual J 3 and J 5 
model deviation signals shown in Figures 4.6 and 4.8. It is plotted with the actual J 3 
model deviation signal for reference. 



Figure 5.33. Effective and Actual / 3 Model Deviation Signals. 
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Thus, the comparisons of interest are between the true effective J 2 and J 2 variations and 
the estimated J 2 and7 3 variations. Continuing with the results for the previous boxcar 
case, Figures 5.34 and 5.35 show the true effective and estimated J 2 and the effective J 2 
residuals respectively. 



Figure 5.34. True Effective J 2 and Estimated J 2 Using Boxcars with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 
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Figure 5.35. Effective J 2 Residuals Using Boxcars with J A and J 5 Model Deviation 
Signals Present but not Estimated. 

Clearly, the boxcar estimates for J 2 agree quite well with the effective values as 
expected. Figures 5.36 and 5.37 show the true effective and estimated / 3 and the 
effective 7 3 residuals respectively. 



Figure 5.36. True Effective J 3 and Estimated J 3 Using Boxcars with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 
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Figure 5.37. Effective 7 3 Residuals Using Boxcars with 7 4 and J 5 Model Deviation 
Signals Present but not Estimated. 

The results for 7 3 also agree quite well with the effective values. Table 5.6 summarizes 
the RMS statistics for the C t , 7 2 , and7 3 residuals. The RMS statistics of the respective 
model deviation signals are also shown for reference. 


Table 5.6. RMS for C„ J 2 , and 7 3 Residuals Using Boxcars with 7 4 and 7 5 Model 
Deviation Signals Present but not Estimated. 


Parameter 

Model Deviation 
Signal RMS 

Residual RMS 
Using Boxcars 

Q 

1.37 x 10' 12 m/s 2 

2.24 x 10" 13 m/s 2 

J 2 

2.85 x lO' 10 

9.24 x 10‘" 

h 

3.10 x lO' 10 

1.03 x 10' 10 


Figures 5.38 to 5.40 show theRTN respective residuals. 







Figure 5.39. Transverse Residuals Using Boxcars with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 




Months Past 1 January 1986 


Figure 5.41. 3-d Position Residuals Using Boxcars with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 


Figures C. 1 through C.6 in Appendix C show the orbit residuals in terms of the 
Keplerian elements (semi-major axis a, eccentricity e, inclination i, argument of 
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periapse CD, longitude of ascending node Q, and argument of latitude co +/, where / is 
true anomaly). 

Finally, the range residuals are determined. Figure 5.42 shows the range 
residuals, and Table 5.7 summarizes the RMS statistics for the RTN, 3-d position, and 
range residuals for the boxcar mode. 


10 
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Months Past 1 January 1986 


Figure 5.42. Range Residuals Using Boxcars with J 4 and J 5 Model Deviation Signals 
Present but not Estimated. 
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Table 5.7. RMS for RTN and Range Residuals Using Boxcars with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 


Residual 

RMS (cm) 

Radial 

0.81 

Transverse 

6.79 

Normal 

3.79 

3-d Position 

7.82 

Range 

2.38 


These positional residuals and range residuals are slightly larger but statistically 
similar to those obtained from the boxcar solution where no J 4 or J s model deviation 
signal is present in the data (refer to Table 5.2). These slightly larger residuals indicate 
that the model deviation signals in the J 4 and J 5 parameters have not been completely 
removed. These residuals, along with the individual C t , J 2 , and J 3 parameter residuals, 
are now compared to the process noise SRIF solution. 

5.3.2 Process Noise SRIF 

The process noise SRIF is implemented with stochastic estimates made for C„ 
J 2 , and J 3 , but not for J 4 or J 5 even though model deviation signals exist for J 4 and J 5 . 
Stochastic estimates of each parameter, based on specific values of r and ex, are made at 
every time that an observation exists. The values of t and o used for each parameter 
are the same as those that were used in the previous process noise case where no J 4 or 
J 5 model deviation signal is present in the data (refer to Table 5.3). 
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The stochastic estimates for the C t , J 2 , and 7 3 temporal variations are compared 
to the true variations. Figure 5.43 shows the true C, and the estimated C, and Figure 
5.44 shows the C t residuals. 



Figure 5.43. True Q and Estimated C, Using Process Noise with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 
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Figure 5.44. C, Residuals Using Process Noise withJ 4 and7 5 Model Deviation Signals 
Present but not Estimated. 


The stochastic estimates of Q appear just as superb as in the previous stochastic case 
where no J 4 or J 5 model deviation signals are present. Figures 5.45 and 5.46 show the 
true effective and estimated J 2 and the effective J 2 residuals respectively. 



Figure 5.45. True Effective J 2 and Estimated J 2 Using Process Noise with J 4 and J s 
Model Deviation Signals Present but not Estimated. 
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Figure 5.46. Effective J 2 Residuals Using Process Noise with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 


Figures 5.46 and 5.48 show the true effective and estimated J 3 and the effective / 3 
residuals respectively. 



Figure 5.47. True Effective / 3 and Estimated J 3 Using Process Noise with J 4 and J 5 
Model Deviation Signals Present but not Estimated. 
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Figure 5.48. Effective J 3 Residuals Using Process Noise with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 


Table 5.9 summarizes the RMS statistics for the Q, J 2 , and J 3 residuals using process 
noise. The respective residuals from the boxcar solution are shown for reference (from 
Table 5.6). The RMS statistics of the respective model deviation signals are also 
shown for reference. 


Table 5.8. RMS for C„ J 2 , and J 3 Residuals Using Process Noise and Boxcars with / 4 
and J 5 Model Deviation Signals Present but not Estimated. 


Parameter 

Model Deviation 
Signal RMS 

Residual RMS 
Using Boxcars 

Residual RMS 
Using Process Noise 

Q 

1.37 x 10' 12 m/s 2 

2.24 x 10' 13 m/s 2 

6.30 x 10~' 4 m/s 2 

h 

2.85 x 10' 10 

9.24 x lO - ' 1 

1.82 x lO’" 

h 

3.10 x 10' 10 

1.03 x 10 -, ° 

2.46 x lO' 11 
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Again, the process noise residuals for C„ J 2 , and J 3 are noticeably improved relative to 
the respective boxcar residuals. 

Next, the RTN residuals are shown in Figures 5.49 to 5.51. The 3-d position 
residuals are shown in Figure 5.52. 


1.0 
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Figure 5.49. Radial Residuals Using Process Noise with / 4 and J 5 Model Deviation 
Signals Present but not Estimated. 
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Figure 5.51. Normal Residuals Using Process Noise with J A and J s Model Deviation 
Signals Present but not Estimated. 
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Figure 5.52. 3-d Position Residuals Using Process Noise with / 4 and J 5 Model 

Deviation Signals Present but not Estimated. 


Figures D.l through D.6 in Appendix D show the orbit residuals in terms of the 
Keplerian elements (semi-major axis a, eccentricity e, inclination i, argument of 

periapse a), longitude of ascending node Q, and argument of latitude (o +/, where / is 

true anomaly). 

Figure 5.53 shows the range residuals, and Table 5.10 summarizes the RMS 
statistics for the RTN, 3-d position, and range residuals (with comparison to the 
respective boxcar statistics from Table 5.7 for reference). 
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Months Past 1 January 1986 

Figure 5.53. Range Residuals Using Process Noise with J 4 and J s Model Deviation 
Signals Present but not Estimated. 


Table 5.9. RMS for RTN and Range Residuals Using Process Noise and Boxcars 
with 7 4 and J 5 Model Deviation Signals Present but not Estimated. 


Residual 

RMS Using 
Boxcars (cm) 

RMS Using 
Process Noise (cm) 

Radial 

0.81 

0.55 

Transverse 

6.79 

0.82 

Normal 

3.79 

0.79 

3-d Position 

7.82 

1.27 

Range 

2.38 

1.12 


Again, relative to the process noise results shown in Table 5.5 where no J 4 or J 5 model 
deviation signal exists in the data, these process noise results show slightly larger 
residuals. These slightly larger residuals indicate that the model deviation signals in the 
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J 4 and J s parameters have not been completely removed. It is also clear that these 
process noise residuals are much better than the corresponding boxcar residuals. 

Figures 5.54 to 5.58 show the RTN, 3-d position, and range residuals from 
both the standard boxcar SRIF and the process noise SRIF solutions side by side on 
the same scale for visual comparison. 



Figure 5.54. Radial Residuals Using Boxcars (a) and Process Noise (b) with J 4 and J 5 
Model Deviation Signals Present but not Estimated. 
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Figure 5.55. Transverse Residuals Using Boxcars (a) and Process Noise (b) with 7 4 
and J 5 Model Deviation Signals Present but not Estimated. 



0369 12 0369 12 

Months Past 1 January 1986 Months Past 1 January 1986 


(a) (b) 

Figure 5.56. Normal Residuals Using Boxcars (a) and Process Noise (b) with J 4 and 
J s Model Deviation Signals Present but not Estimated. 
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Months Past 1 January 1986 


Months Past 1 January 1986 


(a) 


(b) 


Figure 5.57. 3-d Position Residuals Using Boxcars (a) and Process Noise (b) with J 4 
and J s Model Deviation Signals Present but not Estimated. 



Figure 5.58. Range Residuals Using Boxcars (a) and Process Noise (b) with J 4 and J 5 
Model Deviation Signals Present but not Estimated. 
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Again, the stochastic filtering method produces a temporal solution for C„ 
effective J 2 , and effective J 3 which is more accurate than the standard boxcar solution. 
The RMS for the C t , J 2 , and J 3 residuals are reduced to approximately 25% of the 
respective residuals using the standard boxcar method (a 75% improvement). 
Similarly, by using the stochastic filtering method the total 3-d positional residual is 
reduced to approximately 15% of the 3-d positional residual resulting from the standard 
boxcar filtering method (an 85% improvement). Once again, the range residuals are 
reduced by about 50%. The improvement in the orbit accuracy and temporal resolution 
of the model deviation signals is clearly noticeable when using stochastic process noise 
parameters to estimate the model deviation signals. 

5.4 Three Year Arc: C t , J 2 , and J 3 Model Deviation Signals 

A simulated three year data arc with model deviation signals present in C„ J 2 , 
and A and no model deviation signals present in / 4 or J 5 is processed. The three year 
arc is only processed using the stochastic process noise filtering method. A non- 
stochastic solution for the three year data arc requires substantially more time 
computationally due to the significantly larger satellite state, and is not determined in 
this study. Thus, the three year stochastic solution is presented alone without 
comparison to a standard non-stochastic solution. The three year data arc is processed 
with a single estimate made for the satellite state ( r and r ), and no correction made for 
fi, / 4 , or J 5 (since no model deviation signal is present in the ji, i 4 , or J 5 parameters, no 
correction is needed). For this three year process noise SRIF solution, two preliminary 
corrections are made to the satellite state ( r and r ) while fixing all other parameters. 
This gives a reference orbit that is sufficiently close to the true orbit over the entire three 
year arc, thus minimizing linearization errors. A third and final iteration is performed 
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using the process noise mode once these preliminary, non-stochastic corrections are 
made to the initial conditions of the satellite. 

5.4.1 True Model Deviation Signals for the Three Year Arc 

As detailed in section 4.4, specific model deviation signals are introduced into 
the model parameters as the simulated SLR measurements are generated. For the three 
year arc, model deviation signals are introduced into the C„ J 2 , and i 3 parameters. No 
model deviation signal is present in J 4 or J s for the three year arc. The model deviation 
signals for the three year arc are derived in the same manner that those from the one 
year arc are derived (refer to section 4.4). The only difference is that the model 
deviation signals for the three year arc are based on C„ J 2 , and 7 3 estimates from the 
three year period from January 1, 1986 to December 31, 1988. The true model 
deviation signals introduced into the three year simulated data arc are now detailed. 

5.4.1.1 Drag Model Deviation Signal 

For the along track drag parameter C„ the model deviation signal shown in 
Figure 5.59 is added to the nominal C, value of -3.5 picometer/s 2 during the simulation 
of the observations. 
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Figure 5.59. True Model Deviation Signal for the C t Parameter Based on Current Best 
Estimates. 

This model deviation signal is again taken from Tapley et al [1993] which gives 15 day 
estimates for the observed along track drag for LAGEOS over a 14 year period. The 15 
day estimates falling between January 1, 1986 and December 31, 1988 are also shown 
in Figure 5.59 for reference. This three years of interpolated along track drag variation, 
with the average value subtracted out, is what is shown in Figure 5.59. 

5.4.1.2 J 2 and J 3 Model Deviation Signal 

As was the case for the one year arc, the model deviation signal for the J 2 
coefficient consists of a secular and nonsecular term. The model deviation signal for 
the / 3 coefficient is purely nonsecular. The nonsecular part of the J 2 and J 3 model 
deviation signal is again taken from Nerem etal. [1993] which gives monthly estimates 
for the nonsecular variations of J 2 and J 3 over the time period from 1980 to 1989. Their 
monthly estimates for the nonsecular variations from January, 1986 to December, 1988 
are interpolated with a natural cubic spline. This three years of interpolated J 2 and J 3 
nonsecular variation, with the secular variation for J 2 added back, is what is shown in 




J 2 Model Deviation Signal 


Figures 5.60 and 5.61 with the monthly estimates they are based on shown for 
reference (unnormalized). This is the total model deviation signal (total temporal 
variation) added to the nominal JGM-2 J 2 and J 3 values. 



Months Past 1 January 1986 


Figure 5.60. True Model Deviation Signal for the J 2 Parameter Based on Current Best 
Estimates. 
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Figure 5.61. True Model Deviation Signal for the 7 3 Parameter Based on Current Best 
Estimates. 


5.4.2 Process Noise SRIF 

The process noise SRIF is implemented with stochastic estimates made for C„ 
J 2 , and J 3 . Stochastic estimates of each parameter are made at every time that an 
observation exists during the three year arc. The values of t and o used for each 
parameter are the same as those that were used for the one year arc stochastic solutions 
(refer to Table 5.3). 

The stochastic estimates for the C t , J 2 , and J 3 temporal variations are compared 
to the truth. Figure 5.62 shows the true C, and the estimated C, and Figure 5.63 shows 
the C, residuals, or the true C, minus the estimated C,. 
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Figure 5.62. True C t and Estimated C t Using Process Noise. 



Figure 5.63. C t Residuals Using Process Noise. 

The three year stochastic estimate of Q appears excellent. No noticeable degradation in 
the solution is observed relative to the previous stochastic estimate for C t from the one 
year arc. Figures 5.64 and 5.65 show the true J 2 and estimated J 2 and the J 2 residuals 
respectively. 
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Figure 5.64. True J 2 and Estimated J 2 Using Process Noise. 



Figure 5.65. J 2 Residuals Using Process Noise. 

Likewise, Figures 5.66 and 5.67 show the true J 3 and estimated J 3 and the / 3 residuals 
respectively. 
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Figure 5.66. True / 3 and Estimated / 3 Using Process Noise. 



Figure 5.67. / 3 Residuals Using Process Noise. 

The stochastic estimates of J 2 and J 3 are also excellent. Relative to the previous 
stochastic estimates for J 2 and J 3 from the one year arc, the three year stochastic 
estimates appear to have similar temporal resolution and accuracy. Table 5.10 
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summarizes the RMS statistics for the C„ J 2 , and J 3 residuals. The RMS statistics of 
the respective model deviation signals are shown for reference. 

Table 5.10. RMS for C„ J 2 , and7 3 Residuals Using Process Noise. 


Parameter 

Model Deviation 

Residual RMS 


Signal RMS 

Using Process Noise 

c t 

1.10 x 10“ 12 m/s 2 

4.65 x 10 14 m/s 2 

J 2 

2.22 x 10" 10 

1.17 x 10“" 

y 3 

3.14 x lO' 10 

2.06 x 10-” 


The three year process noise residuals for C„ J 2 , and 7 3 are slightly improved relative to 
the respective one year process noise residuals (refer to Table 5.4). 

Next, the RTN positional residuals are shown in Figures 5.68 to 5.70. The 
total RTN, or 3-d position residuals are shown in Figure 5.71. Figures 5.68 to 5.71 
indicate the presence of numerical integration errors at the end of the three year arc. 
This can be dealt with by using integration methods which minimize the accumulation 
of this type of error. 







Figure 5.69. Transverse Residuals Using Process Noise. 






3-d Position Residuals (cm) 



Figure 5.71. 3-d Position Residuals Using Process Noise. 


Figure 5.72 shows the range residuals, and Table 5.11 summarizes the RMS 


statistics for the RTN, 3-d position, and range residuals. Again, since 1 cm RMS noise 
is present in the simulated observations, it is expected that the RMS of the range 
residuals will approach 1 cm as the estimated solution approaches the truth. 
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Months Past 1 January 1986 
Figure 5.72. Range Residuals Using Process Noise. 


Table 5.11. RMS for RTN and Range Residuals Using Process Noise. 


Residual 

RMS Using 
Process Noise (cm) 

Radial 

0.58 

Transverse 

0.49 

Normal 

0.54 

3-d Position 

0.93 

Range 

1.14 


The stochastic filtering method produces a three year temporal solution for C,, 
J 2 , and i 3 which is comparable to the one year stochastic solution. The RMS for the 


stochastic C t , J 2 , and J 3 residuals are slightly less for the three year arc than the one year 
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arc. The orbit position and range residuals are similar for both the three year and one 
year stochastic solutions as well (refer to Table 5.5). Reducing the errors associated 
with numerically integrating the orbit for three years may lead to improved orbit 
position and range residuals for the three year solution. 



CHAPTER 6 


CONCLUSIONS 

6.1 Summary and Discussion 

The feasibility of determining the temporal variations in geodynamically 
interesting parameters, such as the low degree coefficients of the Earth's gravity field, 
through the use of LAGEOS SLR tracking data has been analyzed. A simulation that 
included realistic variations in J 2 and J 3 , and also the LAGEOS along track drag effect, 
was carried out to evaluate the capability of a stochastic filter to track these variations 
using the relatively sparse SLR data. 

Various conclusions can be drawn in assessing the results of this analysis. 
Overall, the filtering results have shown that a stochastic filter can accurately track the 
temporal variations in LAGEOS along track drag, as well as in the J 2 and J 3 gravity 
field coefficients. And the accuracy of these estimates is such that the expected 
variations in these parameters are readily observable. 

In addition, these (stochastic) filtering results are found to provide much better 
accuracy and much better temporal resolution as compared to the conventional (boxcar) 
estimation procedure. These improvements are with respect to results obtained from a 
standard, non-stochastic filter making semi-monthly estimates for C and monthly 
estimates for J 2 and 7 3 . Similarly, the positional accuracy of the orbit is improved by 
using the process noise filter. These improvements are observed both when estimating 
the parameters with and without J 4 and J 5 variations present. Generally, by using the 
process noise filtering approach, the residual RMS errors for the variations in C„ J 2 , 
and J 3 are reduced to approximately 25% of the residual RMS errors obtained using the 
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standard non-stochastic filtering approach. This demonstrates that it is potentially 
feasible to use process noise filtering techniques to improve the temporal resolutions 
and accuracies of LAGEOS along track drag variations and second and third degree 
zonal harmonic coefficient variations. It is also feasible to improve the orbit positional 
accuracy by using multiple stochastic parameters which estimate specific geophysical 
parameters as opposed to simply estimating additional accelerations with these 
parameters. By using the process noise filter, the residual RMS errors for the total 3-d 
orbit position are reduced to approximately 15% of the residual RMS errors obtained 
using the standard filter. Thus, both the orbit positional accuracy and the temporal 
resolution of geodynamic parameters can be improved simultaneously. These 
improvements in the estimates of the C„ J 2 , and J 3 parameters and the orbit itself 
translate into a 50% reduction in the RMS error in the range residuals. 

In addition, the stochastic results from the three year arc show that these 
improvements are possible for longer arcs. While no comparisons to a non-stochastic 
solution for the three year arc are made, the accuracies of the estimates for the model 
parameters and the orbit position are similar to those from the one year stochastic 
solution. This suggests that the ability of the stochastic filter to improve the orbit 
position and temporal resolution of geodynamic parameters simultaneously may be 
possible for longer arcs (such as three years or more) as well. In fact, the accuracies of 
the stochastic estimates for the C„ J 2 , and / 3 parameters were slightly better for the three 
year arc than those from the one year arc. However, other factors such as numerical 
integration errors do become more important for the three year arc, as evidenced by the 
orbit position and range residuals near the end of the three year simulation case. 
Nevertheless, this type of error can be handled effectively, and the accuracies of the 
three year stochastic estimates for the C„ J 2 , and / 3 parameters are very encouraging. 
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It seems reasonable to conclude that improving the temporal resolution of 
similar parameters is possible. Examples would be solar radiation pressure coefficient 
variations, higher degree effective geopotential coefficient variations, and 
nonconservative force coefficients in general. Most variations in Earth orientation 
parameters can be estimated with constants at a daily resolution (daily boxcars). The 
strength in determining these parameters geometrically allows estimates to be made at a 
much higher frequency. However, resolution may still be improved with stochastic 
filtering methods. 

It also seems feasible to extend this stochastic filtering approach to other 
satellites. Satellites whose orbits are reasonably predictable and known and whose 
observations are highly accurate with noise levels which are suitably low would be 
appropriate candidates. If the noise level associated with the observations is too high, it 
may mask the temporal variations of interest. Depending on the particular satellite and 
the impact of specific geodynamic parameters on its dynamics, the feasibility of 
estimating these parameters stochastically with any degree of accuracy may or may not 
be possible. 

6.1.1 Benefits of Stochastic Filtering 

Some general benefits of process noise filtering are discussed as they relate to 
this research and possibly other applications. It is clear that stochastic filters have the 
potential to be used to improve orbit accuracies and temporal resolutions of specific 
geophysical parameters. In addition to improving orbit accuracies along with temporal 
estimates of geophysical parameters, another major benefit is computational in nature. 
In using a standard filter to estimate multiple boxcar estimates over a long arc, 
additional columns of a portion of the state transition matrix must be integrated since the 
state vector requires an additional parameter for each boxcar. As the length of the data 


f 
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arc increases, the number of boxcars required for each parameter increases in order to 
maintain a given temporal resolution for each parameter. For example, in estimating 12 
monthly boxcar estimates for a particular parameter throughout a one year arc rather 
than a single year long estimate adds 1 1 columns to the partition of the state transition 
matrix which corresponds to the satellite accelerations (velocity rows of the state 
vector). For this study, estimating 15 day boxcars for along track drag and monthly 
estimates for J 2 and J 3 over the one year arc requires 48 such columns as opposed to 
three columns if a single year long estimate is made. With the stochastic filter, 
however, no additional columns of the state transition matrix need to be integrated since 
the size of the satellite state vector remains fixed. Only one parameter is required to 
estimate a specific temporal variation over any arc length since independent estimates 
are possible at each observation time. For arcs that are one year or longer, this results 
in a truly significant time savings computationally. Due to this computational 
disadvantage for the standard non-stochastic approach, a boxcar solution for the three 
year arc was not computed for comparison to the stochastic solution. In order to 
estimate a single dynamically consistent orbit for a three year boxcar solution, 72 
parameters for C t , 36 parameters for J 2 , and 36 parameters for / 3 would need to be 
estimated simultaneously. Compared to the stochastic method, the computational cost 
is extraordinary. Moreover, the temporal resolution of the model parameters and 
accuracy of the orbit position is likely to be no better than those from the one year 
boxcar solution. 

Further, for variations which are observable in the presence of the measurement 
noise and temporally resolved equally by both filtering approaches, the stochastic 
approach is simply a more elegant way of estimating those variations that the standard 
approach must estimate with boxcars. By using process noise parameters to estimate 
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unpredictable variations in geodynamic parameters stochastically, no increase in the size 
of the satellite state vector is required. By specifying the time correlation window for 
the variation, controlled by r, and the amplitude of the variation, controlled by cr, the 
need to split up the arc into multiple boxcars is eliminated. While some innovative 
algorithms may minimize the computational hindrances associated with the standard 
boxcar filter to some extent, the process noise approach is inherently computationally 
advantageous with respect to maintaining a smaller state vector and state transition 
matrix while providing higher temporal resolution. In addition, with respect to boxcar 
estimates, a very high temporal resolution may be achieved through the definition of T 
and cr since separate estimates are made at every time an observation is accumulated. 

6.1.2 Drawbacks of Stochastic Filtering 

One drawback related to stochastic filtering observed in this study relates to 
iterating in order to improve the reference trajectory. Improving the initial reference 
trajectory by iterating and producing an improved reference trajectory with stochastic 
solutions of specific state parameters through the use of a process noise filter is 
generally not desirable. If a stochastic solution were to be used for the new reference 
trajectory, then new <r parameters would need to be used for each iteration since 
corrections to the stochastic solution would likely be a different order of magnitude than 
the previous correction. While quite straightforward theoretically, the determination of 
successive values for <7 for each process noise parameter that produce a meaningful 
improved stochastic solution may not be practical. In addition, the use of a stochastic 
solution for a reference trajectory may complicate the filtering depending on the 
robustness of the filter in propagating a reference trajectory and the particular 
application. This approach is definitely possible, but not likely to be advantageous. 
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This minor drawback is quite trivial since determining a suitable reference trajectory 
without using stochastics is generally not a problem. 

Another drawback is the use of stochastics to estimate completely unknown 
variations in a parameter. Not knowing some basic trends regarding the amplitude and 
frequency of the variations’ deviation from an average or reference value is 
unfavorable, particularly if multiple parameters are involved. While estimates can be 
made with no such insight, they may not be reasonable until a proper r and o are 
chosen. It is possible, however, to estimate a parameter with a random walk process 
by simply choosing an appropriate variance q for the parameter (assuming the variation 
is expected to be continuous). In this study, random walk solutions, with accuracies 
comparable to the colored noise solutions previously presented, were generated using a 
variance which was on the order of magnitude of the square of the amplitude of the 
model deviation signal for the parameter. Based on a random walk solution, additional 
insights into a proper rand care likely, since information regarding the amplitude and 
frequency of the unknown model deviation signal is gained. The practical significance 
of this drawback is minor since a preliminary standard boxcar solution for a given 
variation will likely provide adequate insights into the temporal behavior of the 
variation. 

Finally, the possibilities of using stochastic filtering to improve orbit accuracies 
and temporal resolutions of geodynamic parameters should be kept in perspective. 
While this study has shown that it is feasible to make such improvements, it does not 
suggest that applying process noise filtering to all satellite solutions will improve every 
aspect of each solution. Although it is very encouraging to verify via simulated 
LAGEOS SLR data that variations in parameters such as along track drag, J 2 , and / 3 
can be recovered and estimated more accurately with process noise parameters while 
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improving the orbit position accuracy, overly optimistic extensions of these findings 
without similar verification is improper. 

6.2 Recommendations for Future Work 

Some recommendations for future studies are now summarized. A similar 
simulation using a different satellite state may be of interest. By using a satellite state 
composed of Keplerian elements, variations from a reference secularly precessing 
ellipse might be estimated with or instead of the actual geodynamic variations. This 
would lead to direct comparisons of those variations to temporal variations in the 
geopotential or other model variations. Further insights as to the relationships between 
these variations may be gained. 

Another study which assesses the advantages of a hybrid standard-stochastic 
filter may prove valuable. Filters using both boxcars and stochastic parameters to 
estimate particular variations would be more robust. Estimating some parameters, such 
as daily polar motion parameters or tracking station positions and movements, with 
boxcars, and other parameters, such as gravity or non-conservative force variations, 
stochastically would be interesting. Also, by incorporating data from multiple 
satellites, the determination of true variations in higher degree geopotential coefficients 
might be possible since aliasing could be reduced. This multi-satellite study might use 
a simulation to assess the feasibility of determining multiple variations of geopotential 
coefficients. This type of study would contribute to improvements in the modeling of 
the static geopotential field by better accommodating the temporal variations in gravity. 
Since many of the long-wavelength gravity coefficients are only known to the same 
level as their observed temporal variations [Nerem et al., 1993], improvements in the 
estimates of these coefficients must involve dealing with the temporal variations in an 
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appropriate manner. A similar study addressing the estimation of temporal variations in 
the Mars gravity field is also recommended. 

Finally, a study using actual data, perhaps LAGEOS, to compare a stochastic 
filtering solution to a standard boxcar solution over a multi-year arc is recommended. 
Knowing that it is feasible to estimate parameters such as C t , J 2 , and J 3 stochastically, 
such a study would seem productive. The dynamical model would need to be more 
complex, such as that used in actual LAGEOS filtering. The stochastic estimates could 
then be compared to the existing boxcar estimates of the corresponding parameters. 
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APPENDIX A 


KEPLERIAN RESIDUALS FOR ONE YEAR BOXCAR ARC WITH C„ J 2 , AND 7 3 

MODEL DEVIATION SIGNALS 


The following Figures are orbit differences (with respect to the truth) resulting 
from the one year arc with C„ J 2 , and / 3 model deviation signals present using the 
standard boxcar method. They are shown in terms of the Keplerian elements (semi- 
major axis a , eccentricity e, inclination i, argument of periapse a), longitude of 

ascending node Q, and argument of latitude co +/, where /is true anomaly). 
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Figure A. 1 . Semi -Major Axis Residuals Using Boxcars. 
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Figure A.2. Eccentricity Residuals Using Boxcars. 
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Figure A.3. Inclination Residuals Using Boxcars. 
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Figure A.5. Longitude of Ascending Node Residuals Using Boxcars. 
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Figure A.6. Argument of Latitude Residuals Using Boxcars. 




APPENDIX B 


KEPLERIAN RESIDUALS FOR ONE YEAR PROCESS NOISE ARC WITH C„ J 2 , 
AND / 3 MODEL DEVIATION SIGNALS 


The following Figures are orbit differences (with respect to the truth) resulting 
from the one year arc with C„ J 2 , and J 2 model deviation signals present using the 
process noise method. They are shown in terms of the Keplerian elements (semi-major 

axis a, eccentricity e, inclination i, argument of periapse 0 ), longitude of ascending node 


Q, and argument of latitude 0) +/, where /is true anomaly). 



Figure B. 1. Semi -Major Axis Residuals Using Process Noise. 
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Months Past 1 January 1986 


Figure B.2. Eccentricity Residuals Using Process Noise. 



Figure B.3. Inclination Residuals Using Process Noise. 
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k Argument of Periapse Residuals Using Process Noise. 



Figure B.5. Longitude of Ascending Node Residuals Using Process Noise. 






Months Past 1 January 1986 


Figure B.6. Argument of Latitude Residuals Using Process Noise. 





APPENDIX C 


KEPLERIAN RESIDUALS FOR ONE YEAR BOXCAR ARC WITH C„ J 2 , / 3 , J 4 
AND J s MODEL DEVIATION SIGNALS 

The following Figures are orbit differences (with respect to the truth) resulting 
from the one year arc with C„ J 2 , 7 3 , 7 4 , and J s model deviation signals present using 
the standard boxcar method. They are shown in terms of the Keplerian elements (semi- 

major axis a, eccentricity e, inclination i, argument of periapse co, longitude of 
ascending node Q, and argument of latitude (O +/, where/ is true anomaly). 
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Figure C. 1 . Semi -Major Axis Residuals Using Boxcars with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 
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Figure C.2. Eccentricity Residuals Using Boxcars with / 4 and J 5 Model Deviation 
Signals Present but not Estimated. 
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Figure C.3. Inclination Residuals Using Boxcars with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 
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Figure C.4. Argument of Periapse Residuals Using Boxcars with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 



Figure C.5. Longitude of Ascending Node Residuals Using Boxcars with J 4 and J 5 
Model Deviation Signals Present but not Estimated. 
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Figure C.6. Argument of Latitude Residuals Using Boxcars with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 





APPENDIX D 


KEPLERIAN RESIDUALS FOR ONE YEAR PROCESS NOISE ARC WITH C t , J 2 , 
7 3 , y 4 AND J 5 MODEL DEVIATION SIGNALS 


The following Figures are orbit differences (with respect to the truth) resulting 
from the one year arc with C„ J 2 , J 2 , i 4 , and J 5 model deviation signals present using 
the process noise method. They are shown in terms of the Keplerian elements (semi- 
major axis a, eccentricity e, inclination i, argument of periapse (o, longitude of 


ascending node Q, and argument of latitude a> +/, where /is true anomaly). 
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Figure D.l. Semi-Major Axis Residuals Using Process Noise with J 4 and J 5 Model 
Deviation Signals Present but not Estimated. 
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Figure D.2. Eccentricity Residuals Using Process Noise with J 4 and J s Model 
Deviation Signals Present but not Estimated. 
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Figure D.3. Inclination Residuals Using Process Noise with J 4 and J 5 Model Deviation 
Signals Present but not Estimated. 
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Months Past 1 January 1986 


Argument of Periapse Residuals Using Process Noise with J A and J s 
Nation Signals Present but not Estimated. 
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Figure D.5. Longitude of Ascending Node Residuals Using Process Noise with J 4 and 
J 5 Model Deviation Signals Present but not Estimated. 
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Figure D.6. Argument of Latitude Residuals Using Process Noise with J 4 and J 5 
Model DeviationSignals Present but not Estimated. 



